From: Benjamin Auder <benjamin.auder@somewhere> Date: Thu, 23 Feb 2017 17:47:41 +0000 (+0100) Subject: Drop useless files (we will use R 'parallel' package) X-Git-Url: https://git.auder.net/variants/img/current/assets/css/pieces/cp.svg?a=commitdiff_plain;h=7f1a6cf08a4d4d67e8a95b8c1c0cc74ff3deb5a4;p=valse.git Drop useless files (we will use R 'parallel' package) --- diff --git a/src/test/generate_test_data/helpers/constructionModelesLassoMLE.R b/R/constructionModelesLassoMLE.R similarity index 100% rename from src/test/generate_test_data/helpers/constructionModelesLassoMLE.R rename to R/constructionModelesLassoMLE.R diff --git a/src/test/generate_test_data/helpers/constructionModelesLassoRank.R b/R/constructionModelesLassoRank.R similarity index 100% rename from src/test/generate_test_data/helpers/constructionModelesLassoRank.R rename to R/constructionModelesLassoRank.R diff --git a/src/adapters/a.constructionModelesLassoMLE.c b/src/adapters/a.constructionModelesLassoMLE.c deleted file mode 100644 index ec519a9..0000000 --- a/src/adapters/a.constructionModelesLassoMLE.c +++ /dev/null @@ -1,92 +0,0 @@ -#include <R.h> -#include <Rdefines.h> -#include "constructionModelesLassoMLE.h" - -SEXP constructionModelesLassoMLE( - SEXP phiInit_, - SEXP rhoInit_, - SEXP piInit_, - SEXP gamInit_, - SEXP mini_, - SEXP maxi_, - SEXP gamma_, - SEXP glambda_, - SEXP X_, - SEXP Y_, - SEXP seuil_, - SEXP tau_, - SEXP A1_, - SEXP A2_ -) { - // Get matrices dimensions - int n = INTEGER(getAttrib(X_, R_DimSymbol))[0]; - SEXP dim = getAttrib(phiInit_, R_DimSymbol); - int p = INTEGER(dim)[0]; - int m = INTEGER(dim)[1]; - int k = INTEGER(dim)[2]; - int L = length(glambda_); - - //////////// - // INPUTS // - //////////// - - // get scalar parameters - int mini = INTEGER_VALUE(mini_); - int maxi = INTEGER_VALUE(maxi_); - double gamma = NUMERIC_VALUE(gamma_); - double seuil = NUMERIC_VALUE(seuil_); - double tau = NUMERIC_VALUE(tau_); - - // Get pointers from SEXP arrays ; WARNING: by columns ! - double* phiInit = REAL(phiInit_); - double* rhoInit = REAL(rhoInit_); - double* piInit = REAL(piInit_); - double* gamInit = REAL(gamInit_); - double* glambda = REAL(glambda_); - double* X = REAL(X_); - double* Y = REAL(Y_); - int* A1 = INTEGER(A1_); - int* A2 = INTEGER(A2_); - - ///////////// - // OUTPUTS // - ///////////// - - 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; - PROTECT(dimRho = allocVector(INTSXP, 4)); - int* pDimRho = INTEGER(dimRho); - pDimRho[0] = m; pDimRho[1] = m; pDimRho[2] = k; pDimRho[3] = L; - PROTECT(phi = allocArray(REALSXP, dimPhi)); - PROTECT(rho = allocArray(REALSXP, dimRho)); - PROTECT(pi = allocMatrix(REALSXP, k, L)); - PROTECT(llh = allocMatrix(REALSXP, L, 2)); - double *pPhi=REAL(phi), *pRho=REAL(rho), *pPi=REAL(pi), *pllh=REAL(llh); - - ///////////////////////////////////////// - // Call to constructionModelesLassoMLE // - ///////////////////////////////////////// - - constructionModelesLassoMLE_core( - phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2, - 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", "llh"}; //lists labels - PROTECT(listNames = allocVector(STRSXP,4)); - for (int i=0; i<4; 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, rho); - SET_VECTOR_ELT(listParams, 2, pi); - SET_VECTOR_ELT(listParams, 3, llh); - - UNPROTECT(8); - return listParams; -} diff --git a/src/adapters/a.constructionModelesLassoRank.c b/src/adapters/a.constructionModelesLassoRank.c deleted file mode 100644 index 0e069d4..0000000 --- a/src/adapters/a.constructionModelesLassoRank.c +++ /dev/null @@ -1,79 +0,0 @@ -#include <R.h> -#include <Rdefines.h> -#include "constructionModelesLassoRank.h" - -SEXP constructionModelesLassoRank( - SEXP Pi_, - SEXP Rho_, - SEXP mini_, - SEXP maxi_, - SEXP X_, - SEXP Y_, - SEXP tau_, - SEXP A1_, - SEXP rangmin_, - SEXP rangmax_ -) { - // Get matrices dimensions - SEXP dimX = getAttrib(X_, R_DimSymbol); - int n = INTEGER(dimX)[0]; - int p = INTEGER(dimX)[1]; - SEXP dimRho = getAttrib(Rho_, R_DimSymbol); - int m = INTEGER(dimRho)[0]; - int k = INTEGER(dimRho)[2]; - int L = INTEGER(getAttrib(A1_, R_DimSymbol))[1]; - - //////////// - // INPUTS // - //////////// - - // get scalar parameters - int mini = INTEGER_VALUE(mini_); - int maxi = INTEGER_VALUE(maxi_); - double tau = NUMERIC_VALUE(tau_); - double rangmin = NUMERIC_VALUE(rangmin_); - double rangmax = NUMERIC_VALUE(rangmax_); - - // Get pointers from SEXP arrays ; WARNING: by columns ! - double* Pi = REAL(Pi_); - double* Rho = REAL(Rho_); - double* X = REAL(X_); - double* Y = REAL(Y_); - int* A1 = INTEGER(A1_); - - ///////////// - // OUTPUTS // - ///////////// - - int Size = pow(rangmax-rangmin+1,k); - 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(llh = allocMatrix(REALSXP, L*Size, 2)); - double *pPhi=REAL(phi), *pllh=REAL(llh); - - ////////////////////////////////////////// - // Call to constructionModelesLassoRank // - ////////////////////////////////////////// - - constructionModelesLassoRank_core( - Pi,Rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax, - 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", "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, llh); - - UNPROTECT(5); - return listParams; -} diff --git a/src/adapters/a.selectiontotale.c b/src/adapters/a.selectiontotale.c deleted file mode 100644 index 3bfab9f..0000000 --- a/src/adapters/a.selectiontotale.c +++ /dev/null @@ -1,90 +0,0 @@ -#include <R.h> -#include <Rdefines.h> -#include "selectiontotale.h" - -SEXP selectiontotale( - SEXP phiInit_, - SEXP rhoInit_, - SEXP piInit_, - SEXP gamInit_, - SEXP mini_, - SEXP maxi_, - SEXP gamma_, - SEXP glambda_, - SEXP X_, - SEXP Y_, - SEXP seuil_, - SEXP tau_ -) { - // Get matrices dimensions - SEXP dimX = getAttrib(X_, R_DimSymbol); - int n = INTEGER(dimX)[0]; - int p = INTEGER(dimX)[1]; - SEXP dimRhoInit = getAttrib(rhoInit_, R_DimSymbol); - int m = INTEGER(dimRhoInit)[0]; - int k = INTEGER(dimRhoInit)[2]; - int L = length(glambda_); - - //////////// - // INPUTS // - //////////// - - // get scalar parameters - int mini = INTEGER_VALUE(mini_); - int maxi = INTEGER_VALUE(maxi_); - double gamma = NUMERIC_VALUE(gamma_); - double seuil = NUMERIC_VALUE(seuil_); - double tau = NUMERIC_VALUE(tau_); - - // Get pointers from SEXP arrays ; WARNING: by columns ! - double* phiInit = REAL(phiInit_); - double* rhoInit = REAL(rhoInit_); - double* piInit = REAL(piInit_); - double* gamInit = REAL(gamInit_); - double* glambda = REAL(glambda_); - double* X = REAL(X_); - double* Y = REAL(Y_); - - ///////////// - // OUTPUTS // - ///////////// - - 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(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)); - int *pA1=INTEGER(A1), *pA2=INTEGER(A2); - double *pRho=REAL(rho), *pPi=REAL(pi); - - ///////////////////////////// - // Call to selectiontotale // - ///////////////////////////// - - selectiontotale_core( - phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau, - pA1,pA2,pRho,pPi, - n,p,m,k,L); - - // Build list from OUT params and return it - SEXP listParams, listNames; - PROTECT(listParams = allocVector(VECSXP, 4)); - char* lnames[4] = { "A1", "A2", "rho", "pi" }; //lists labels - PROTECT(listNames = allocVector(STRSXP, 4)); - for (int i=0; i<4; i++) - SET_STRING_ELT(listNames,i,mkChar(lnames[i])); - setAttrib(listParams, R_NamesSymbol, listNames); - 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/constructionModelesLassoMLE.c b/src/sources/constructionModelesLassoMLE.c deleted file mode 100644 index 34e5808..0000000 --- a/src/sources/constructionModelesLassoMLE.c +++ /dev/null @@ -1,214 +0,0 @@ -#include "EMGLLF.h" -#include "utils.h" -#include <stdlib.h> -#include <gsl/gsl_linalg.h> -#include <omp.h> - -// TODO: comment on constructionModelesLassoMLE purpose -void constructionModelesLassoMLE_core( - // IN parameters - const Real* phiInit, // parametre initial de moyenne renormalisé - const Real* rhoInit, // parametre initial de variance renormalisé - const Real* piInit,// parametre initial des proportions - const Real* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon - int mini,// nombre minimal d'itérations dans l'algorithme EM - int maxi,// nombre maximal d'itérations dans l'algorithme EM - Real gamma,// valeur de gamma : puissance des proportions dans la pénalisation - //pour un Lasso adaptatif - const Real* glambda, // valeur des paramètres de régularisation du Lasso - const Real* X, // régresseurs - const Real* Y, // réponse - Real seuil,// seuil pour prendre en compte une variable - Real tau,// seuil pour accepter la convergence - const int* A1, // matrice des coefficients des parametres selectionnes - const int* A2, // matrice des coefficients des parametres non selectionnes - // OUT parameters - Real* phi,// estimateur ainsi calculé par le Lasso - Real* rho,// estimateur ainsi calculé par le Lasso - Real* pi, // estimateur ainsi calculé par le Lasso - Real* llh, // estimateur ainsi calculé par le Lasso - // additional size parameters - int n, // taille de l'echantillon - int p, // nombre de covariables - int m, // taille de Y (multivarié) - int k, // nombre de composantes - int L) // taille de glambda -{ - //preparation: phi,rho,pi = 0, llh=+Inf - for (int u=0; u<p*m*k*L; u++) - phi[u] = 0.; - for (int u=0; u<m*m*k*L; u++) - rho[u] = 0.; - for (int u=0; u<k*L; u++) - pi[u] = 0.; - for (int u=0; u<L*2; u++) - llh[u] = INFINITY; - - //initiate parallel section - int lambdaIndex; - omp_set_num_threads(OMP_NUM_THREADS); - #pragma omp parallel default(shared) private(lambdaIndex) - { - #pragma omp for schedule(dynamic,CHUNK_SIZE) nowait - for (lambdaIndex=0; lambdaIndex<L; lambdaIndex++) - { - //a = A1[,1,lambdaIndex] ; a = a[a!=0] - int* a = (int*)malloc(p*sizeof(int)); - int lengthA = 0; - for (int j=0; j<p; j++) - { - if (A1[ai(j,0,lambdaIndex,p,m+1,L)] != 0) - a[lengthA++] = A1[ai(j,0,lambdaIndex,p,m+1,L)] - 1; - } - if (lengthA == 0) - { - free(a); - continue; - } - - //Xa = X[,a] - Real* Xa = (Real*)malloc(n*lengthA*sizeof(Real)); - for (int i=0; i<n; i++) - { - for (int j=0; j<lengthA; j++) - Xa[mi(i,j,n,lengthA)] = X[mi(i,a[j],n,p)]; - } - - //phia = phiInit[a,,] - Real* phia = (Real*)malloc(lengthA*m*k*sizeof(Real)); - for (int j=0; j<lengthA; j++) - { - for (int mm=0; mm<m; mm++) - { - for (int r=0; r<k; r++) - phia[ai(j,mm,r,lengthA,m,k)] = phiInit[ai(a[j],mm,r,p,m,k)]; - } - } - - //Call to EMGLLF - Real* phiLambda = (Real*)malloc(lengthA*m*k*sizeof(Real)); - Real* rhoLambda = (Real*)malloc(m*m*k*sizeof(Real)); - Real* piLambda = (Real*)malloc(k*sizeof(Real)); - Real* LLF = (Real*)malloc((maxi+1)*sizeof(Real)); - Real* S = (Real*)malloc(lengthA*m*k*sizeof(Real)); - EMGLLF_core(phia,rhoInit,piInit,gamInit,mini,maxi,gamma,0.,Xa,Y,tau, - phiLambda,rhoLambda,piLambda,LLF,S, - n,lengthA,m,k); - free(Xa); - free(phia); - free(LLF); - free(S); - - //Assign results to current variables - for (int j=0; j<lengthA; j++) - { - for (int mm=0; mm<m; mm++) - { - for (int r=0; r<k; r++) - phi[ai4(a[j],mm,r,lambdaIndex,p,m,k,L)] = phiLambda[ai(j,mm,r,lengthA,m,k)]; - } - } - free(phiLambda); - for (int u=0; u<m; u++) - { - for (int v=0; v<m; v++) - { - for (int r=0; r<k; r++) - rho[ai4(u,v,r,lambdaIndex,m,m,k,L)] = rhoLambda[ai(u,v,r,m,m,k)]; - } - } - free(rhoLambda); - for (int r=0; r<k; r++) - pi[mi(r,lambdaIndex,k,L)] = piLambda[r]; - free(piLambda); - - int dimension = 0; - int* b = (int*)malloc(m*sizeof(int)); - for (int j=0; j<p; j++) - { - //b = A2[j,2:dim(A2)[2],lambdaIndex] ; b = b[b!=0] - int lengthB = 0; - for (int mm=0; mm<m; mm++) - { - if (A2[ai(j,mm+1,lambdaIndex,p,m+1,L)] != 0) - b[lengthB++] = A2[ai(j,mm+1,lambdaIndex,p,m+1,L)] - 1; - } - if (lengthB > 0) - { - //phi[A2[j,1,lambdaIndex],b,,lambdaIndex] = 0. - for (int mm=0; mm<lengthB; mm++) - { - for (int r=0; r<k; r++) - phi[ai4(A2[ai(j,0,lambdaIndex,p,m+1,L)]-1, b[mm], r, lambdaIndex, p, m, k, L)] = 0.; - } - } - - //c = A1[j,2:dim(A1)[2],lambdaIndex] ; dimension = dimension + sum(c!=0) - for (int mm=0; mm<m; mm++) - { - if (A1[ai(j,mm+1,lambdaIndex,p,m+1,L)] != 0) - dimension++; - } - } - free(b); - - int signum; - Real* densite = (Real*)calloc(L*n,sizeof(Real)); - Real sumLogDensit = 0.0; - gsl_matrix* matrix = gsl_matrix_alloc(m, m); - gsl_permutation* permutation = gsl_permutation_alloc(m); - Real* YiRhoR = (Real*)malloc(m*sizeof(Real)); - Real* XiPhiR = (Real*)malloc(m*sizeof(Real)); - for (int i=0; i<n; i++) - { - for (int r=0; r<k; r++) - { - //compute det(rho(:,:,r,lambdaIndex)) [TODO: avoid re-computations] - for (int u=0; u<m; u++) - { - for (int v=0; v<m; v++) - matrix->data[u*m+v] = rho[ai4(u,v,r,lambdaIndex,m,m,k,L)]; - } - gsl_linalg_LU_decomp(matrix, permutation, &signum); - Real detRhoR = gsl_linalg_LU_det(matrix, signum); - - //compute Y(i,:)*rho(:,:,r,lambdaIndex) - for (int u=0; u<m; u++) - { - YiRhoR[u] = 0.0; - for (int v=0; v<m; v++) - YiRhoR[u] += Y[mi(i,v,n,m)] * rho[ai4(v,u,r,lambdaIndex,m,m,k,L)]; - } - - //compute X(i,a)*phi(a,:,r,lambdaIndex) - for (int u=0; u<m; u++) - { - XiPhiR[u] = 0.0; - for (int v=0; v<lengthA; v++) - XiPhiR[u] += X[mi(i,a[v],n,p)] * phi[ai4(a[v],u,r,lambdaIndex,p,m,k,L)]; - } - // NOTE: On peut remplacer X par Xa dans ce dernier calcul, - // mais je ne sais pas si c'est intéressant ... - - // compute dotProduct < delta . delta > - Real dotProduct = 0.0; - for (int u=0; u<m; u++) - dotProduct += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]); - - densite[mi(lambdaIndex,i,L,n)] += - (pi[mi(r,lambdaIndex,k,L)]*detRhoR/pow(sqrt(2.0*M_PI),m))*exp(-dotProduct/2.0); - } - sumLogDensit += log(densite[lambdaIndex*n+i]); - } - llh[mi(lambdaIndex,0,L,2)] = sumLogDensit; - llh[mi(lambdaIndex,1,L,2)] = (dimension+m+1)*k-1; - - free(a); - free(YiRhoR); - free(XiPhiR); - free(densite); - gsl_matrix_free(matrix); - gsl_permutation_free(permutation); - } - } -} diff --git a/src/sources/constructionModelesLassoMLE.h b/src/sources/constructionModelesLassoMLE.h deleted file mode 100644 index 06e500d..0000000 --- a/src/sources/constructionModelesLassoMLE.h +++ /dev/null @@ -1,34 +0,0 @@ -#ifndef valse_constructionModelesLassoMLE_H -#define valse_constructionModelesLassoMLE_H - -#include "utils.h" - -void constructionModelesLassoMLE_core( - // IN parameters - const Real* phiInit, - const Real* rhoInit, - const Real* piInit, - const Real* gamInit, - int mini, - int maxi, - Real gamma, - const Real* glambda, - const Real* X, - const Real* Y, - Real seuil, - Real tau, - const int* A1, - const int* A2, - // OUT parameters - Real* phi, - Real* rho, - Real* pi, - Real* llh, - // additional size parameters - int n, - int p, - int m, - int k, - int L); - -#endif diff --git a/src/sources/constructionModelesLassoRank.c b/src/sources/constructionModelesLassoRank.c deleted file mode 100644 index 1115311..0000000 --- a/src/sources/constructionModelesLassoRank.c +++ /dev/null @@ -1,138 +0,0 @@ -#include <stdlib.h> -#include <omp.h> -#include <gsl/gsl_linalg.h> -#include "EMGrank.h" -#include "utils.h" - -// TODO: comment on constructionModelesLassoRank purpose -void constructionModelesLassoRank_core( - // IN parameters - const Real* pi,// parametre initial des proportions - const Real* rho, // parametre initial de variance renormalisé - int mini, // nombre minimal d'itérations dans l'algorithme EM - int maxi, // nombre maximal d'itérations dans l'algorithme EM - const Real* X,// régresseurs - const Real* Y,// réponse - Real tau, // seuil pour accepter la convergence - const int* A1, // matrice des coefficients des parametres selectionnes - int rangmin, //rang minimum autorisé - int rangmax, //rang maximum autorisé - // OUT parameters (all pointers, to be modified) - Real* phi,// estimateur ainsi calculé par le Lasso - Real* llh,// estimateur ainsi calculé par le Lasso - // additional size parameters - int n,// taille de l'echantillon - int p,// nombre de covariables - int m,// taille de Y (multivarié) - int k,// nombre de composantes - int L)// taille de glambda -{ - //On cherche les rangs possiblement intéressants - int deltaRank = rangmax-rangmin+1; - int Size = (int)pow(deltaRank,k); - int* Rank = (int*)malloc(Size*k*sizeof(int)); - for (int r=0; r<k; r++) - { - // On veut le tableau de toutes les combinaisons de rangs possibles - // Dans la première colonne : on répète (rangmax-rangmin)^(k-1) chaque chiffre : - // ça remplit la colonne - // Dans la deuxieme : on répète (rangmax-rangmin)^(k-2) chaque chiffre, - // et on fait ça (rangmax-rangmin)^2 fois - // ... - // Dans la dernière, on répète chaque chiffre une fois, - // et on fait ça (rangmin-rangmax)^(k-1) fois. - int indexInRank = 0; - int value = 0; - while (indexInRank < Size) - { - int upperBoundU = (int)pow(deltaRank,k-r-1); - for (int u=0; u<upperBoundU; u++) - Rank[mi(indexInRank++,r,Size,k)] = rangmin + value; - value = (value+1) % deltaRank; - } - } - - //Initialize phi to zero, because unactive variables won't be assigned - for (int i=0; i<p*m*k*L*Size; i++) - phi[i] = 0.0; - for (int i=0; i<L*Size*2; i++) - llh[i] = INFINITY; - - //initiate parallel section - int lambdaIndex; - omp_set_num_threads(OMP_NUM_THREADS); - #pragma omp parallel default(shared) private(lambdaIndex) - { - #pragma omp for schedule(dynamic,CHUNK_SIZE) nowait - for (lambdaIndex=0; lambdaIndex<L; lambdaIndex++) - { - //On ne garde que les colonnes actives : active sera l'ensemble des variables informatives - int* active = (int*)malloc(p*sizeof(int)); - int longueurActive = 0; - for (int j=0; j<p; j++) - { - if (A1[mi(j,lambdaIndex,p,L)] != 0) - active[longueurActive++] = A1[mi(j,lambdaIndex,p,L)] - 1; - } - if (longueurActive == 0) - continue; - - //from now on, longueurActive > 0 - Real* phiLambda = (Real*)malloc(longueurActive*m*k*sizeof(Real)); - Real LLF; - for (int j=0; j<Size; j++) - { - int* rank = (int*)malloc(k*sizeof(int)); - for (int r=0; r<k; r++) - rank[r] = Rank[mi(j,r,Size,k)]; - Real* Xactive = (Real*)malloc(n*longueurActive*sizeof(Real)); - for (int i=0; i<n; i++) - { - for (int jj=0; jj<longueurActive; jj++) - Xactive[mi(i,jj,n,longueurActive)] = X[mi(i,active[jj],n,p)]; - } - Real* piLambda = (Real*)malloc(k*sizeof(Real)); - for (int r=0; r<k; r++) - piLambda[r] = pi[mi(r,lambdaIndex,k,L)]; - Real* rhoLambda = (Real*)malloc(m*m*k*sizeof(Real)); - for (int u=0; u<m; u++) - { - for (int v=0; v<m; v++) - { - for (int r=0; r<k; r++) - rhoLambda[ai(u,v,r,m,m,k)] = rho[ai4(u,v,r,lambdaIndex,m,m,k,L)]; - } - } - EMGrank_core(piLambda,rhoLambda,mini,maxi,Xactive,Y,tau,rank, - phiLambda,&LLF, - n,longueurActive,m,k); - free(rank); - free(Xactive); - free(piLambda); - free(rhoLambda); - //llh[(lambdaIndex-1)*Size+j,] = c(LLF, ...) - llh[mi(lambdaIndex*Size+j,0,L*Size,2)] = LLF; - //sum(Rank[j,] * (length(active)- Rank[j,] + m)) ) - Real dotProduct = 0.0; - for (int r=0; r<k; r++) - dotProduct += Rank[mi(j,r,Size,k)] * (longueurActive-Rank[mi(j,r,Size,k)]+m); - llh[mi(lambdaIndex*Size+j,1,Size*L,2)] = dotProduct; - //phi[active,,,(lambdaIndex-1)*Size+j] = res$phi - for (int jj=0; jj<longueurActive; jj++) - { - for (int mm=0; mm<m; mm++) - { - for (int r=0; r<k; r++) - { - phi[ai4(active[jj],mm,r,lambdaIndex*Size+j,p,m,k,L*Size)] = - phiLambda[ai(jj,mm,r,longueurActive,m,k)]; - } - } - } - } - free(active); - free(phiLambda); - } - } - free(Rank); -} diff --git a/src/sources/constructionModelesLassoRank.h b/src/sources/constructionModelesLassoRank.h deleted file mode 100644 index dbe8f59..0000000 --- a/src/sources/constructionModelesLassoRank.h +++ /dev/null @@ -1,29 +0,0 @@ -#ifndef valse_constructionModelesLassoRank_H -#define valse_constructionModelesLassoRank_H - -#include "utils.h" - -// Main job on raw inputs (after transformation from mxArray) -void constructionModelesLassoRank_core( - // IN parameters - const Real* Pi, - const Real* Rho, - int mini, - int maxi, - const Real* X, - const Real* Y, - Real tau, - const int* A1, - int rangmin, - int rangmax, - // OUT parameters - Real* phi, - Real* llh, - // additional size parameters - int n, - int p, - int m, - int k, - int L); - -#endif diff --git a/src/sources/selectiontotale.c b/src/sources/selectiontotale.c deleted file mode 100644 index 8c1b768..0000000 --- a/src/sources/selectiontotale.c +++ /dev/null @@ -1,143 +0,0 @@ -#include <stdlib.h> -#include <omp.h> -#include <math.h> -#include "EMGLLF.h" -#include "utils.h" - -// Main job on raw inputs (after transformation from mxArray) -void selectiontotale_core( - // IN parameters - const Real* phiInit, // parametre initial de moyenne renormalisé - const Real* rhoInit, // parametre initial de variance renormalisé - const Real* piInit,// parametre initial des proportions - const Real* gamInit, // paramètre initial des probabilités a posteriori de chaque échantillon - int mini, // nombre minimal d'itérations dans lambdaIndex'algorithme EM - int maxi, // nombre maximal d'itérations dans lambdaIndex'algorithme EM - Real gamma, // valeur de gamma : puissance des proportions dans la pénalisation pour un Lasso adaptatif - const Real* glambda, // valeur des paramètres de régularisation du Lasso - const Real* X,// régresseurs - const Real* Y,// réponse - Real seuil, // seuil pour prendre en compte une variable - Real tau, // seuil pour accepter la convergence - // OUT parameters (all pointers, to be modified) - int* A1, // matrice des coefficients des parametres selectionnes - int* A2, // matrice des coefficients des parametres non selectionnes - Real* Rho,// estimateur ainsi calculé par le Lasso - Real* Pi,// estimateur ainsi calculé par le Lasso - // additional size parameters - int n,// taille de lambdaIndex'echantillon - int p,// nombre de covariables - int m,// taille de Y (multivarié) - int k,// nombre de composantes - int L) // taille de glambda -{ - // Fill outputs with zeros: they might not be assigned - for (int u=0; u<p*(m+1)*L; u++) - { - A1[u] = 0; - A2[u] = 0; - } - for (int u=0; u<m*m*k*L; u++) - Rho[u] = 0.0; - for (int u=0; u<k*L; u++) - Pi[u] = 0.0; - - //initiate parallel section - int lambdaIndex; - omp_set_num_threads(OMP_NUM_THREADS); - #pragma omp parallel default(shared) private(lambdaIndex) - { - #pragma omp for schedule(dynamic,CHUNK_SIZE) nowait - for (lambdaIndex=0; lambdaIndex<L; lambdaIndex++) - { - //allocate output variables - 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)); - EMGLLF_core(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda[lambdaIndex],X,Y,tau, - phi,rho,pi,LLF,S, - n,p,m,k); - free(LLF); - free(S); - - //Si un des coefficients est supérieur au seuil, on garde cette variable - int* selectedVariables = (int*)calloc(p*m,sizeof(int)); - int* discardedVariables = (int*)calloc(p*m,sizeof(int)); - int atLeastOneSelectedVariable = 0; - for (int j=0; j<p; j++) - { - int cpt = 0; - int cpt2 = 0; - for (int jj=0; jj<m; jj++) - { - Real maxPhi = 0.0; - for (int r=0; r<k; r++) - { - if (fabs(phi[ai(j,jj,r,p,m,k)]) > maxPhi) - maxPhi = fabs(phi[ai(j,jj,r,p,m,k)]); - } - if (maxPhi > seuil) - { - selectedVariables[mi(j,cpt,p,m)] = jj+1; - atLeastOneSelectedVariable = 1; - cpt++; - } - else - { - discardedVariables[mi(j,cpt2,p,m)] = jj+1; - cpt2++; - } - } - } - free(phi); - - if (atLeastOneSelectedVariable) - { - int* vec = (int*)malloc(p*sizeof(int)); - int vecSize = 0; - for (int j=0; j<p; j++) - { - if (selectedVariables[mi(j,0,p,m)] != 0) - vec[vecSize++] = j; - } - - //A1 - for (int j=0; j<p; j++) - A1[ai(j,0,lambdaIndex,p,m+1,L)] = (j < vecSize ? vec[j]+1 : 0); - for (int j=0; j<vecSize; j++) - { - for (int jj=1; jj<=m; jj++) - A1[ai(j,jj,lambdaIndex,p,m+1,L)] = selectedVariables[mi(vec[j],jj-1,p,m)]; - } - //A2 - for (int j=0; j<p; j++) - A2[ai(j,0,lambdaIndex,p,m+1,L)] = j+1; - for (int j=0; j<p; j++) - { - for (int jj=1; jj<=m; jj++) - A2[ai(j,jj,lambdaIndex,p,m+1,L)] = discardedVariables[mi(j,jj-1,p,m)]; - } - //Rho - for (int j=0; j<m; j++) - { - for (int jj=0; jj<m; jj++) - { - for (int r=0; r<k; r++) - Rho[ai4(j,jj,r,lambdaIndex,m,m,k,L)] = rho[ai(j,jj,r,m,m,k)]; - } - } - //Pi - for (int r=0; r<k; r++) - Pi[mi(r,lambdaIndex,k,L)] = pi[r]; - free(vec); - } - - free(selectedVariables); - free(discardedVariables); - free(rho); - free(pi); - } - } -} diff --git a/src/sources/selectiontotale.h b/src/sources/selectiontotale.h deleted file mode 100644 index fea6e39..0000000 --- a/src/sources/selectiontotale.h +++ /dev/null @@ -1,33 +0,0 @@ -#ifndef valse_selectiontotale_H -#define valse_selectiontotale_H - -#include "utils.h" - -// Main job on raw inputs (after transformation from mxArray) -void selectiontotale_core( - // IN parameters - const Real* phiInit, - const Real* rhoInit, - const Real* piInit, - const Real* gamInit, - int mini, - int maxi, - Real gamma, - const Real* glambda, - const Real* X, - const Real* Y, - Real seuil, - Real tau, - // OUT parameters - int* A1, - int* A2, - Real* Rho, - Real* Pi, - // additional size parameters - int n, - int p, - int m, - int k, - int L); - -#endif diff --git a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.R b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.R deleted file mode 100644 index dbce359..0000000 --- a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoMLE.R +++ /dev/null @@ -1,65 +0,0 @@ -source("helpers/constructionModelesLassoMLE.R") - -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)) -{ - tau = 1e-6; - seuil = 1e-15; - L = length(glambda); - A1 = array(0, dim=c(p, m+1, L)) - A2 = array(0, dim=c(p, m+1, L)) - for (i in 1:L) - { - A1[,1,i] = seq_len(p) - A1[1:5,2,i] = seq_len(5) - A2[,1,i] = seq_len(p) - A2[1:5,2,i] = seq_len(5) - } - require(valse) - params = valse:::basicInitParameters(n, p, m, k) - xy = valse:::generateXYdefault(n, p, m, k) - - testFolder = "../data/" - dir.create(testFolder, showWarnings=FALSE, mode="0755") - #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(glambda), paste(testFolder,"glambda",sep=""), - row.names=F, col.names=F) - write.table(as.double(xy$X), paste(testFolder,"X",sep=""), - 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(seuil), paste(testFolder,"seuil",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.double(A2), paste(testFolder,"A2",sep=""), - row.names=F, col.names=F) - write.table(as.integer(c(n,p,m,k,L)), paste(testFolder,"dimensions",sep=""), - row.names=F, col.names=F) - - res = constructionModelesLassoMLE( - params$phiInit,params$rhoInit,params$piInit,params$gamInit, - mini,maxi,gamma,glambda,xy$X,xy$Y,seuil,tau,A1,A2) - - #save output - 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_constructionModelesLassoRank.R b/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoRank.R deleted file mode 100644 index 559bd9e..0000000 --- a/src/test/generate_test_data/generateRunSaveTest_constructionModelesLassoRank.R +++ /dev/null @@ -1,49 +0,0 @@ -source("helpers/constructionModelesLassoRank.R") - -generateRunSaveTest_constructionModelesLassoRank = function(n=200, p=15, m=10, k=3, L=12, mini=5, - maxi=10, gamma=1.0, rangmin=3, rangmax=6) -{ - tau = 1e-6 - pi = matrix(1./k, nrow=k, ncol=L) - rho = array(dim=c(m,m,k,L)) - for (l in 1:L) - { - for (r in 1:k) - rho[,,r,l] = diag(1,m) - } - A1 = matrix(seq_len(p), nrow=p, ncol=L) - require(valse) - xy = valse:::generateXYdefault(n, p, m, k) - - testFolder = "../data/" - dir.create(testFolder, showWarnings=FALSE, mode="0755") - #save inputs - write.table(as.double(pi), paste(testFolder,"pi",sep=""), - row.names=F, col.names=F) - write.table(as.double(rho),paste(testFolder,"rho",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(xy$X), paste(testFolder,"X",sep=""), - 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=""), - 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,L)),paste(testFolder,"dimensions",sep=""), - row.names=F, col.names=F) - - res = constructionModelesLassoRank(pi,rho,mini,maxi,xy$X,xy$Y,tau,A1,rangmin,rangmax) - - #save output - 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_selectiontotale.R b/src/test/generate_test_data/generateRunSaveTest_selectiontotale.R deleted file mode 100644 index 0ea03bf..0000000 --- a/src/test/generate_test_data/generateRunSaveTest_selectiontotale.R +++ /dev/null @@ -1,49 +0,0 @@ -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)) -{ - tau = 1e-6; - seuil = 1e-15; - require(valse) - params = valse:::basicInitParameters(n, p, m, k) - xy = valse:::generateXYdefault(n, p, m, k) - L = length(glambda) - - testFolder = "data/" - dir.create(testFolder, showWarnings=FALSE, mode="0755") - #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(glambda), paste(testFolder,"lambda",sep=""), - row.names=F, col.names=F) - write.table(as.double(xy$X), paste(testFolder,"X",sep=""), - 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(seuil), paste(testFolder,"seuil",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,L)), paste(testFolder,"dimensions",sep=""), - row.names=F, col.names=F) - - res = selectiontotale(params$phiInit,params$rhoInit,params$piInit,params$gamInit, - mini,maxi,gamma,glambda,xy$X,xy$Y,seuil, tau) - - #save output - 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/test.constructionModelesLassoMLE.c b/src/test/test.constructionModelesLassoMLE.c deleted file mode 100644 index b607467..0000000 --- a/src/test/test.constructionModelesLassoMLE.c +++ /dev/null @@ -1,83 +0,0 @@ -#include "constructionModelesLassoMLE.h" -#include "test_utils.h" -#include <stdlib.h> - -#include <stdio.h> - -int main(int argc, char** argv) -{ - 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); - - //////////// - // 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"); - //////////// - - ///////////// - // 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)); - ///////////// - - ///////////////////////////////////////// - // Call to constructionModelesLassoMLE // - constructionModelesLassoMLE_core( - phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,glambda,X,Y,seuil,tau,A1,A2, - phi,rho,pi,llh, - n,p,m,k,L); - ///////////////////////////////////////// - - free(phiInit); - free(rhoInit); - free(piInit); - free(gamInit); - free(X); - free(Y); - free(A1); - free(A2); - free(glambda); - - // Compare to reference outputs - Real* ref_phi = readArray_real("phi"); - compareArray_real("phi", phi, ref_phi, p*m*k*L); - free(phi); - free(ref_phi); - - Real* ref_rho = readArray_real("rho"); - compareArray_real("rho", rho, ref_rho, m*m*k*L); - free(rho); - free(ref_rho); - - Real* ref_pi = readArray_real("pi"); - compareArray_real("pi", pi, ref_pi, k*L); - free(pi); - free(ref_pi); - - Real* ref_llh = readArray_real("llh"); - compareArray_real("llh", llh, ref_llh, L*2); - free(llh); - free(ref_llh); - - return 0; -} diff --git a/src/test/test.constructionModelesLassoRank.c b/src/test/test.constructionModelesLassoRank.c deleted file mode 100644 index f60ffea..0000000 --- a/src/test/test.constructionModelesLassoRank.c +++ /dev/null @@ -1,63 +0,0 @@ -#include "constructionModelesLassoRank.h" -#include "test_utils.h" -#include <stdlib.h> -#include <math.h> - -int main(int argc, char** argv) -{ - 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); - - //////////// - // 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"); - //////////// - - ///////////// - // OUTPUTS // - int Size = (int)pow(rangmax-rangmin+1, k); - Real* phi = (Real*)malloc(p*m*k*L*Size*sizeof(Real)); - Real* llh = (Real*)malloc(L*Size*2*sizeof(Real)); - ///////////// - - ///////////////////////////////////////// - // Call to constructionModelesLassoMLE // - constructionModelesLassoRank_core( - pi,rho,mini,maxi,X,Y,tau,A1,rangmin,rangmax, - phi,llh, - n,p,m,k,L); - ///////////////////////////////////////// - - free(rho); - free(pi); - free(X); - free(Y); - free(A1); - - // Compare to reference outputs - Real* ref_phi = readArray_real("phi"); - compareArray_real("phi", phi, ref_phi, p*m*k*L*Size); - free(phi); - free(ref_phi); - - 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 deleted file mode 100644 index 843a04d..0000000 --- a/src/test/test.selectiontotale.c +++ /dev/null @@ -1,77 +0,0 @@ -#include "selectiontotale.h" -#include "test_utils.h" -#include <stdlib.h> - -int main(int argc, char** argv) -{ - 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); - - //////////// - // 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"); - //////////// - - ///////////// - // 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)); - ///////////// - - ///////////////////////////////////////// - // 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(phiInit); - free(rhoInit); - free(piInit); - free(gamInit); - free(glambda); - free(X); - free(Y); - - // Compare to reference outputs - int* ref_A1 = readArray_int("A1"); - compareArray_int("A1", A1, ref_A1, p*(m+1)*L); - free(A1); - free(ref_A1); - - int* ref_A2 = readArray_int("A2"); - compareArray_int("A2", A2, ref_A2, p*(m+1)*L); - free(A2); - free(ref_A2); - - Real* ref_Rho = readArray_real("Rho"); - compareArray_real("Rho", Rho, ref_Rho, m*m*k*L); - free(Rho); - free(ref_Rho); - - Real* ref_Pi = readArray_real("Pi"); - compareArray_real("Pi", Pi, ref_Pi, k*L); - free(Pi); - free(ref_Pi); - - return 0; -}