+++ /dev/null
-#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;
-}
+++ /dev/null
-#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;
-}
+++ /dev/null
-#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;
-}
+++ /dev/null
-#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);
- }
- }
-}
+++ /dev/null
-#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
+++ /dev/null
-#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);
-}
+++ /dev/null
-#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
+++ /dev/null
-#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);
- }
- }
-}
+++ /dev/null
-#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
+++ /dev/null
-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)
-}
+++ /dev/null
-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)
-}
+++ /dev/null
-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)
-}
+++ /dev/null
-#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;
-}
+++ /dev/null
-#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;
-}
+++ /dev/null
-#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;
-}
projects / valse.git / commitdiff
