From 4cab944a7bf905e811d740694c39f16b2b23d7e1 Mon Sep 17 00:00:00 2001 From: Benjamin Auder <benjamin.auder@somewhere> Date: Fri, 2 Dec 2016 02:41:42 +0100 Subject: [PATCH] finished a.EMGLLF.c [to test] ; fix indices in EMGLLF.c [by cols] --- src/Makevars | 2 +- src/adapters/a.EMGLLF.c | 231 +++++----------- src/sources/EMGLLF.c | 374 +++++++++++++------------- src/sources/utils/io.c | 168 ------------ src/sources/utils/io.h | 65 +---- src/sources/utils/omp_num_threads.h | 7 - src/sources/utils/tune_parallelisms.h | 10 + 7 files changed, 266 insertions(+), 591 deletions(-) delete mode 100644 src/sources/utils/io.c delete mode 100644 src/sources/utils/omp_num_threads.h create mode 100644 src/sources/utils/tune_parallelisms.h diff --git a/src/Makevars b/src/Makevars index aa08234..bec17f0 100644 --- a/src/Makevars +++ b/src/Makevars @@ -2,6 +2,6 @@ PKG_CFLAGS=-g -I. PKG_LIBS=-lm -SOURCES = $(wildcard adapters/*.c sources/*.c sources/utils/*.c) +SOURCES = $(wildcard adapters/*.c sources/*.c) OBJECTS = $(SOURCES:.c=.o) diff --git a/src/adapters/a.EMGLLF.c b/src/adapters/a.EMGLLF.c index 241a09f..747a3c9 100644 --- a/src/adapters/a.EMGLLF.c +++ b/src/adapters/a.EMGLLF.c @@ -1,197 +1,86 @@ #include <R.h> #include <Rdefines.h> #include "sources/EMGLLF.h" -#include "sources/utils/io.h" SEXP EMGLLF( - SEXP M_, - SEXP NIix_, - SEXP alpha_, - SEXP h_, - SEXP epsilon_, - SEXP maxiter_, - SEXP symmNeighbs_, - SEXP trace_ + SEXP phiInit_, + SEXP rhoInit_, + SEXP piInit_, + SEXP gamInit_, + SEXP mini_, + SEXP maxi_, + SEXP gamma_, + SEXP lambda_, + SEXP X_, + SEXP Y_, + SEXP tau_ ) { - // get parameters - double alpha = NUMERIC_VALUE(alpha_); - double h = NUMERIC_VALUE(h_); - double epsilon = NUMERIC_VALUE(epsilon_); - int maxiter = INTEGER_VALUE(maxiter_); - int symmNeighbs = LOGICAL_VALUE(symmNeighbs_); - int trace = LOGICAL_VALUE(trace_); - - // extract infos from M and get associate pointer - SEXP dim = getAttrib(M_, R_DimSymbol); - int nrow = INTEGER(dim)[0]; - int ncol = INTEGER(dim)[1]; - // M is always given by columns: easier to process in rows - double* pM = transpose(REAL(M_), nrow, ncol); - - // extract NIix list vectors in a jagged array - int* lengthNIix = (int*)malloc(nrow*sizeof(int)); - int** NIix = (int**)malloc(nrow*sizeof(int*)); - for (int i=0; i<nrow; i++) - { - lengthNIix[i] = LENGTH(VECTOR_ELT(NIix_,i)); - SEXP tmp; - PROTECT(tmp = AS_INTEGER(VECTOR_ELT(NIix_,i))); - NIix[i] = (int*)malloc(lengthNIix[i]*sizeof(int)); - for (int j=0; j<lengthNIix[i]; j++) - NIix[i][j] = INTEGER(tmp)[j]; - UNPROTECT(1); - // WARNING: R indices start at 1, - // so we must lower every index right now to avoid future bug - for (int j=0; j<lengthNIix[i]; j++) - NIix[i][j]--; - } - - // Main call to core algorithm - Parameters params = getVarsWithConvexOptim_core( - pM, lengthNIix, NIix, nrow, ncol, alpha, h, epsilon, maxiter, symmNeighbs, trace); - - // free neighborhoods parameters arrays - free(lengthNIix); - for (int i=0; i<nrow; i++) - free(NIix[i]); - free(NIix); - - // copy matrix F into pF for output to R (1D matrices) - SEXP f; - PROTECT(f = allocMatrix(REALSXP, nrow, ncol)); - double* pF = REAL(f); - for (int i=0; i<nrow; i++) - { - for (int j=0; j<ncol; j++) - pF[i+nrow*j] = params.f[i][j]; - } - // copy theta into pTheta for output to R - SEXP theta; - PROTECT(theta = allocVector(REALSXP, nrow)); - double* pTheta = REAL(theta); - for (int i=0; i<nrow; i++) - pTheta[i] = params.theta[i]; - - // free params.f and params.theta - free(params.theta); - for (int i=0; i<nrow; i++) - free(params.f[i]); - free(params.f); - - // build return list with f and theta - SEXP listParams, listNames; - PROTECT(listParams = allocVector(VECSXP, 2)); - char* lnames[2] = {"f", "theta"}; //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, f); - SET_VECTOR_ELT(listParams, 1, theta); - - UNPROTECT(4); - return listParams; - - - - - - - - // Get matrices dimensions - const mwSize n = mxGetDimensions(prhs[8])[0]; - const mwSize p = mxGetDimensions(prhs[0])[0]; - const mwSize m = mxGetDimensions(prhs[0])[1]; - const mwSize k = mxGetDimensions(prhs[0])[2]; + 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]; //////////// // INPUTS // //////////// - // phiInit - const mwSize* dimPhiInit = mxGetDimensions(prhs[0]); - Real* brPhiInit = matlabToBrArray_real(mxGetPr(prhs[0]), dimPhiInit, 3); - - // rhoInit - const mwSize* dimRhoInit = mxGetDimensions(prhs[1]); - Real* brRhoInit = matlabToBrArray_real(mxGetPr(prhs[1]), dimRhoInit, 3); - - // piInit - Real* piInit = mxGetPr(prhs[2]); - - // gamInit - const mwSize* dimGamInit = mxGetDimensions(prhs[3]); - Real* brGamInit = matlabToBrArray_real(mxGetPr(prhs[3]), dimGamInit, 2); - - // min number of iterations - Int mini = ((Int*)mxGetData(prhs[4]))[0]; - - // max number of iterations - Int maxi = ((Int*)mxGetData(prhs[5]))[0]; - - // gamma - Real gamma = mxGetScalar(prhs[6]); - - // lambda - Real lambda = mxGetScalar(prhs[7]); - - // X - const mwSize* dimX = mxGetDimensions(prhs[8]); - Real* brX = matlabToBrArray_real(mxGetPr(prhs[8]), dimX, 2); - - // Y - const mwSize* dimY = mxGetDimensions(prhs[9]); - Real* brY = matlabToBrArray_real(mxGetPr(prhs[9]), dimY, 2); - - // tau - Real tau = mxGetScalar(prhs[10]); + // get scalar parameters + int mini = INTEGER_VALUE(mini_); + int maxi = INTEGER_VALUE(maxi_); + double gamma = NUMERIC_VALUE(gamma_); + double lambda = NUMERIC_VALUE(lambda_); + 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* X = REAL(X_); + double* Y = REAL(Y_); ///////////// // OUTPUTS // ///////////// - // phi - const mwSize dimPhi[] = {dimPhiInit[0], dimPhiInit[1], dimPhiInit[2]}; - plhs[0] = mxCreateNumericArray(3,dimPhi,mxDOUBLE_CLASS,mxREAL); - Real* phi = mxGetPr(plhs[0]); - - // rho - const mwSize dimRho[] = {dimRhoInit[0], dimRhoInit[1], dimRhoInit[2]}; - plhs[1] = mxCreateNumericArray(3,dimRho,mxDOUBLE_CLASS,mxREAL); - Real* rho = mxGetPr(plhs[1]); - - // pi - plhs[2] = mxCreateNumericMatrix(k,1,mxDOUBLE_CLASS,mxREAL); - Real* pi = mxGetPr(plhs[2]); - - // LLF - plhs[3] = mxCreateNumericMatrix(maxi,1,mxDOUBLE_CLASS,mxREAL); - Real* LLF = mxGetPr(plhs[3]); - - // S - const mwSize dimS[] = {p, m, k}; - plhs[4] = mxCreateNumericArray(3,dimS,mxDOUBLE_CLASS,mxREAL); - Real* S = mxGetPr(plhs[4]); + SEXP phi, rho, pi, LLF, S, dimPhiS, dimRho; + PROTECT(dimPhiS = allocVector(INTSXP, 3)); + int* pDimPhiS = INTEGER(dimPhiS); + pDimPhiS[0] = p; pDimPhiS[1] = m; pDimPhiS[2] = k; + PROTECT(dimRho = allocVector(INTSXP, 3)); + int* pDimRho = INTEGER(dimRho); + pDimRho[0] = m; pDimRho[1] = m; pDimRho[2] = k; + PROTECT(phi = allocArray(REALSXP, dimPhiS)); + PROTECT(rho = allocArray(REALSXP, dimRho)); + PROTECT(pi = allocVector(REALSXP, k)); + PROTECT(LLF = allocVector(REALSXP, maxi-mini+1)); + PROTECT(S = allocArray(REALSXP, dimPhiS)); + double* pPhi=REAL(phi), pRho=REAL(rho), pPi=REAL(pi), pLLF=REAL(LLF), pS=REAL(S); //////////////////// // Call to EMGLLF // //////////////////// - EMGLLF(brPhiInit,brRhoInit,piInit,brGamInit,mini,maxi,gamma,lambda,brX,brY,tau, - phi,rho,pi,LLF,S, + EMGLLF(phiInit,rhoInit,piInit,gamInit,mini,maxi,gamma,lambda,X,Y,tau, + pPhi,pRho,pPi,pLLF,pS, n,p,m,k); - - free(brPhiInit); - free(brRhoInit); - free(brGamInit); - free(brX); - free(brY); - - - - - + // Build list from OUT params and return it + SEXP listParams, listNames; + PROTECT(listParams = allocVector(VECSXP, 5)); + char* lnames[5] = {"phi", "rho", "pi", "LLF", "S"}; //lists labels + PROTECT(listNames = allocVector(STRSXP,5)); + 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, 3, LLF); + SET_ARRAY_ELT(listParams, 4, S); + UNPROTECT(9); + return listParams; } diff --git a/src/sources/EMGLLF.c b/src/sources/EMGLLF.c index 96b81b3..0c39d6e 100644 --- a/src/sources/EMGLLF.c +++ b/src/sources/EMGLLF.c @@ -1,31 +1,31 @@ #include "EMGLLF.h" #include <gsl/gsl_linalg.h> -// TODO: comment on EMGLLF purpose +// TODO: don't recompute indexes every time...... void EMGLLF( // 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 - Real lambda, // valeur du paramètre de régularisation du Lasso - const Real* X, // régresseurs - const Real* Y, // réponse - Real tau, // seuil pour accepter la convergence + const double* phiInit, // parametre initial de moyenne renormalisé + const double* rhoInit, // parametre initial de variance renormalisé + const double* piInit, // parametre initial des proportions + const double* 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 + double gamma, // valeur de gamma : puissance des proportions dans la pénalisation pour un Lasso adaptatif + double lambda, // valeur du paramètre de régularisation du Lasso + const double* X, // régresseurs + const double* Y, // réponse + double tau, // seuil pour accepter la convergence // OUT parameters (all pointers, to be modified) - Real* phi, // parametre de moyenne renormalisé, calculé par l'EM - Real* rho, // parametre de variance renormalisé, calculé par l'EM - Real* pi, // parametre des proportions renormalisé, calculé par l'EM - Real* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres - Real* S, + double* phi, // parametre de moyenne renormalisé, calculé par l'EM + double* rho, // parametre de variance renormalisé, calculé par l'EM + double* pi, // parametre des proportions renormalisé, calculé par l'EM + double* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres + double* S, // additional size parameters - mwSize n, // nombre d'echantillons - mwSize p, // nombre de covariables - mwSize m, // taille de Y (multivarié) - mwSize k) // nombre de composantes dans le mélange + int n, // nombre d'echantillons + int p, // nombre de covariables + int m, // taille de Y (multivarié) + int k) // nombre de composantes dans le mélange { //Initialize outputs copyArray(phiInit, phi, p*m*k); @@ -33,77 +33,77 @@ void EMGLLF( copyArray(piInit, pi, k); zeroArray(LLF, maxi); //S is already allocated, and doesn't need to be 'zeroed' - + //Other local variables //NOTE: variables order is always [maxi],n,p,m,k - Real* gam = (Real*)malloc(n*k*sizeof(Real)); + double* gam = (double*)malloc(n*k*sizeof(double)); copyArray(gamInit, gam, n*k); - Real* b = (Real*)malloc(k*sizeof(Real)); - 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* gam2 = (Real*)malloc(k*sizeof(Real)); - Real* pi2 = (Real*)malloc(k*sizeof(Real)); - Real* Gram2 = (Real*)malloc(p*p*k*sizeof(Real)); - Real* ps = (Real*)malloc(m*k*sizeof(Real)); - Real* nY2 = (Real*)malloc(m*k*sizeof(Real)); - Real* ps1 = (Real*)malloc(n*m*k*sizeof(Real)); - Real* ps2 = (Real*)malloc(p*m*k*sizeof(Real)); - Real* nY21 = (Real*)malloc(n*m*k*sizeof(Real)); - Real* Gam = (Real*)malloc(n*k*sizeof(Real)); - Real* X2 = (Real*)malloc(n*p*k*sizeof(Real)); - Real* Y2 = (Real*)malloc(n*m*k*sizeof(Real)); + double* b = (double*)malloc(k*sizeof(double)); + double* Phi = (double*)malloc(p*m*k*sizeof(double)); + double* Rho = (double*)malloc(m*m*k*sizeof(double)); + double* Pi = (double*)malloc(k*sizeof(double)); + double* gam2 = (double*)malloc(k*sizeof(double)); + double* pi2 = (double*)malloc(k*sizeof(double)); + double* Gram2 = (double*)malloc(p*p*k*sizeof(double)); + double* ps = (double*)malloc(m*k*sizeof(double)); + double* nY2 = (double*)malloc(m*k*sizeof(double)); + double* ps1 = (double*)malloc(n*m*k*sizeof(double)); + double* ps2 = (double*)malloc(p*m*k*sizeof(double)); + double* nY21 = (double*)malloc(n*m*k*sizeof(double)); + double* Gam = (double*)malloc(n*k*sizeof(double)); + double* X2 = (double*)malloc(n*p*k*sizeof(double)); + double* Y2 = (double*)malloc(n*m*k*sizeof(double)); 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)); - Real dist = 0.0; - Real dist2 = 0.0; - Int ite = 0; - Real EPS = 1e-15; - Real* dotProducts = (Real*)malloc(k*sizeof(Real)); - + double* YiRhoR = (double*)malloc(m*sizeof(double)); + double* XiPhiR = (double*)malloc(m*sizeof(double)); + double dist = 0.; + double dist2 = 0.; + int ite = 0; + double EPS = 1e-15; + double* dotProducts = (double*)malloc(k*sizeof(double)); + while (ite < mini || (ite < maxi && (dist >= tau || dist2 >= sqrt(tau)))) { copyArray(phi, Phi, p*m*k); copyArray(rho, Rho, m*m*k); copyArray(pi, Pi, k); - - // Calculs associes a Y et X - for (mwSize r=0; r<k; r++) + + // Calculs associés a Y et X + for (int r=0; r<k; r++) { - for (mwSize mm=0; mm<m; mm++) + for (int mm=0; mm<m; mm++) { //Y2(:,mm,r)=sqrt(gam(:,r)).*transpose(Y(mm,:)); - for (mwSize u=0; u<n; u++) - Y2[u*m*k+mm*k+r] = sqrt(gam[u*k+r]) * Y[u*m+mm]; + for (int u=0; u<n; u++) + Y2[ai(u,mm,r,n,m,k)] = sqrt(gam[mi(u,r,n,k)]) * Y[mi(u,mm,m,n)]; } - for (mwSize i=0; i<n; i++) + for (int i=0; i<n; i++) { //X2(i,:,r)=X(i,:).*sqrt(gam(i,r)); - for (mwSize u=0; u<p; u++) - X2[i*p*k+u*k+r] = sqrt(gam[i*k+r]) * X[i*p+u]; + for (int u=0; u<p; u++) + X2[ai(i,u,r,n,m,k)] = sqrt(gam[mi(i,r,n,k)]) * X[mi(i,u,n,p)]; } - for (mwSize mm=0; mm<m; mm++) + for (int mm=0; mm<m; mm++) { //ps2(:,mm,r)=transpose(X2(:,:,r))*Y2(:,mm,r); - for (mwSize u=0; u<p; u++) + for (int u=0; u<p; u++) { - Real dotProduct = 0.0; - for (mwSize v=0; v<n; v++) - dotProduct += X2[v*p*k+u*k+r] * Y2[v*m*k+mm*k+r]; - ps2[u*m*k+mm*k+r] = dotProduct; + double dotProduct = 0.; + for (int v=0; v<n; v++) + dotProduct += X2[ai(v,u,r,n,m,k)] * Y2[ai(v,mm,r,n,m,k)]; + ps2[ai(u,mm,r,n,m,k)] = dotProduct; } } - for (mwSize j=0; j<p; j++) + for (int j=0; j<p; j++) { - for (mwSize s=0; s<p; s++) + for (int s=0; s<p; s++) { //Gram2(j,s,r)=transpose(X2(:,j,r))*(X2(:,s,r)); - Real dotProduct = 0.0; - for (mwSize u=0; u<n; u++) - dotProduct += X2[u*p*k+j*k+r] * X2[u*p*k+s*k+r]; - Gram2[j*p*k+s*k+r] = dotProduct; + double dotProduct = 0.; + for (int u=0; u<n; u++) + dotProduct += X2[ai(u,j,r,n,p,k)] * X2[ai(u,s,r,n,p,k)]; + Gram2[ai(j,s,r,p,p,k)] = dotProduct; } } } @@ -111,46 +111,46 @@ void EMGLLF( ///////////// // Etape M // ///////////// - + // Pour pi - for (mwSize r=0; r<k; r++) + for (int r=0; r<k; r++) { //b(r) = sum(sum(abs(phi(:,:,r)))); - Real sumAbsPhi = 0.0; - for (mwSize u=0; u<p; u++) - for (mwSize v=0; v<m; v++) - sumAbsPhi += fabs(phi[u*m*k+v*k+r]); + double sumAbsPhi = 0.; + for (int u=0; u<p; u++) + for (int v=0; v<m; v++) + sumAbsPhi += fabs(phi[ai(u,v,r,p,m,k)]); b[r] = sumAbsPhi; } //gam2 = sum(gam,1); - for (mwSize u=0; u<k; u++) + for (int u=0; u<k; u++) { - Real sumOnColumn = 0.0; - for (mwSize v=0; v<n; v++) - sumOnColumn += gam[v*k+u]; + double sumOnColumn = 0.; + for (int v=0; v<n; v++) + sumOnColumn += gam[mi(v,u,n,k)]; gam2[u] = sumOnColumn; } //a=sum(gam*transpose(log(pi))); - Real a = 0.0; - for (mwSize u=0; u<n; u++) + double a = 0.; + for (int u=0; u<n; u++) { - Real dotProduct = 0.0; - for (mwSize v=0; v<k; v++) - dotProduct += gam[u*k+v] * log(pi[v]); + double dotProduct = 0.; + for (int v=0; v<k; v++) + dotProduct += gam[mi(u,v,n,k)] * log(pi[v]); a += dotProduct; } - + //tant que les proportions sont negatives - mwSize kk = 0; + int kk = 0; int pi2AllPositive = 0; - Real invN = 1.0/n; + double invN = 1./n; while (!pi2AllPositive) { //pi2(:)=pi(:)+0.1^kk*(1/n*gam2(:)-pi(:)); - for (mwSize r=0; r<k; r++) + for (int r=0; r<k; r++) pi2[r] = pi[r] + pow(0.1,kk) * (invN*gam2[r] - pi[r]); pi2AllPositive = 1; - for (mwSize r=0; r<k; r++) + for (int r=0; r<k; r++) { if (pi2[r] < 0) { @@ -160,113 +160,113 @@ void EMGLLF( } kk++; } - + //t(m) la plus grande valeur dans la grille O.1^k tel que ce soit décroissante ou constante //(pi.^gamma)*b - Real piPowGammaDotB = 0.0; - for (mwSize v=0; v<k; v++) + double piPowGammaDotB = 0.; + for (int v=0; v<k; v++) piPowGammaDotB += pow(pi[v],gamma) * b[v]; //(pi2.^gamma)*b - Real pi2PowGammaDotB = 0.0; - for (mwSize v=0; v<k; v++) + double pi2PowGammaDotB = 0.; + for (int v=0; v<k; v++) pi2PowGammaDotB += pow(pi2[v],gamma) * b[v]; //transpose(gam2)*log(pi2) - Real prodGam2logPi2 = 0.0; - for (mwSize v=0; v<k; v++) + double prodGam2logPi2 = 0.; + for (int v=0; v<k; v++) prodGam2logPi2 += gam2[v] * log(pi2[v]); while (-invN*a + lambda*piPowGammaDotB < -invN*prodGam2logPi2 + lambda*pi2PowGammaDotB && kk<1000) { //pi2=pi+0.1^kk*(1/n*gam2-pi); - for (mwSize v=0; v<k; v++) + for (int v=0; v<k; v++) pi2[v] = pi[v] + pow(0.1,kk) * (invN*gam2[v] - pi[v]); //pi2 was updated, so we recompute pi2PowGammaDotB and prodGam2logPi2 - pi2PowGammaDotB = 0.0; - for (mwSize v=0; v<k; v++) + pi2PowGammaDotB = 0.; + for (int v=0; v<k; v++) pi2PowGammaDotB += pow(pi2[v],gamma) * b[v]; - prodGam2logPi2 = 0.0; - for (mwSize v=0; v<k; v++) + prodGam2logPi2 = 0.; + for (int v=0; v<k; v++) prodGam2logPi2 += gam2[v] * log(pi2[v]); kk++; } - Real t = pow(0.1,kk); + double t = pow(0.1,kk); //sum(pi+t*(pi2-pi)) - Real sumPiPlusTbyDiff = 0.0; - for (mwSize v=0; v<k; v++) + double sumPiPlusTbyDiff = 0.; + for (int v=0; v<k; v++) sumPiPlusTbyDiff += (pi[v] + t*(pi2[v] - pi[v])); //pi=(pi+t*(pi2-pi))/sum(pi+t*(pi2-pi)); - for (mwSize v=0; v<k; v++) + for (int v=0; v<k; v++) pi[v] = (pi[v] + t*(pi2[v] - pi[v])) / sumPiPlusTbyDiff; - + //Pour phi et rho - for (mwSize r=0; r<k; r++) + for (int r=0; r<k; r++) { - for (mwSize mm=0; mm<m; mm++) + for (int mm=0; mm<m; mm++) { - for (mwSize i=0; i<n; i++) + for (int i=0; i<n; i++) { //< X2(i,:,r) , phi(:,mm,r) > - Real dotProduct = 0.0; - for (mwSize u=0; u<p; u++) - dotProduct += X2[i*p*k+u*k+r] * phi[u*m*k+mm*k+r]; + double dotProduct = 0.0; + for (int u=0; u<p; u++) + dotProduct += X2[ai(i,u,r,n,p,k)] * phi[ai(u,mm,r,n,m,k)]; //ps1(i,mm,r)=Y2(i,mm,r)*dot(X2(i,:,r),phi(:,mm,r)); - ps1[i*m*k+mm*k+r] = Y2[i*m*k+mm*k+r] * dotProduct; - nY21[i*m*k+mm*k+r] = Y2[i*m*k+mm*k+r] * Y2[i*m*k+mm*k+r]; + ps1[ai(i,mm,r,n,m,k)] = Y2[ai(i,mm,r,n,m,k)] * dotProduct; + nY21[ai(i,mm,r,n,m,k)] = Y2[ai(i,mm,r,n,m,k)] * Y2[ai(i,mm,r,n,m,k)]; } //ps(mm,r)=sum(ps1(:,mm,r)); - Real sumPs1 = 0.0; - for (mwSize u=0; u<n; u++) - sumPs1 += ps1[u*m*k+mm*k+r]; - ps[mm*k+r] = sumPs1; + double sumPs1 = 0.0; + for (int u=0; u<n; u++) + sumPs1 += ps1[ai(u,mm,r,n,m,k)]; + ps[mi(mm,r,m,k)] = sumPs1; //nY2(mm,r)=sum(nY21(:,mm,r)); - Real sumNy21 = 0.0; - for (mwSize u=0; u<n; u++) - sumNy21 += nY21[u*m*k+mm*k+r]; - nY2[mm*k+r] = sumNy21; + double sumNy21 = 0.0; + for (int u=0; u<n; u++) + sumNy21 += nY21[ai(u,mm,r,n,m,k)]; + nY2[mi(mm,r,m,k)] = sumNy21; //rho(mm,mm,r)=((ps(mm,r)+sqrt(ps(mm,r)^2+4*nY2(mm,r)*(gam2(r))))/(2*nY2(mm,r))); - rho[mm*m*k+mm*k+r] = ( ps[mm*k+r] + sqrt( ps[mm*k+r]*ps[mm*k+r] - + 4*nY2[mm*k+r] * (gam2[r]) ) ) / (2*nY2[mm*k+r]); + rho[ai(mm,mm,k,m,m,k)] = ( ps[mi(mm,r,m,k)] + sqrt( ps[mi(mm,r,m,k)]*ps[mi(mm,r,m,k)] + + 4*nY2[mi(mm,r,m,k)] * (gam2[r]) ) ) / (2*nY2[mi(mm,r,m,k)]); } } - for (mwSize r=0; r<k; r++) + for (int r=0; r<k; r++) { - for (mwSize j=0; j<p; j++) + for (int j=0; j<p; j++) { - for (mwSize mm=0; mm<m; mm++) + for (int mm=0; mm<m; mm++) { //sum(phi(1:j-1,mm,r).*transpose(Gram2(j,1:j-1,r)))+sum(phi(j+1:p,mm,r).*transpose(Gram2(j,j+1:p,r))) - Real dotPhiGram2 = 0.0; - for (mwSize u=0; u<j; u++) - dotPhiGram2 += phi[u*m*k+mm*k+r] * Gram2[j*p*k+u*k+r]; - for (mwSize u=j+1; u<p; u++) - dotPhiGram2 += phi[u*m*k+mm*k+r] * Gram2[j*p*k+u*k+r]; + double dotPhiGram2 = 0.0; + for (int u=0; u<j; u++) + dotPhiGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)]; + for (int u=j+1; u<p; u++) + dotPhiGram2 += phi[ai(u,mm,r,p,m,k)] * Gram2[ai(j,u,r,p,p,k)]; //S(j,r,mm)=-rho(mm,mm,r)*ps2(j,mm,r)+sum(phi(1:j-1,mm,r).*transpose(Gram2(j,1:j-1,r))) // +sum(phi(j+1:p,mm,r).*transpose(Gram2(j,j+1:p,r))); - S[j*m*k+mm*k+r] = -rho[mm*m*k+mm*k+r] * ps2[j*m*k+mm*k+r] + dotPhiGram2; - if (fabs(S[j*m*k+mm*k+r]) <= n*lambda*pow(pi[r],gamma)) - phi[j*m*k+mm*k+r] = 0; - else if (S[j*m*k+mm*k+r] > n*lambda*pow(pi[r],gamma)) - phi[j*m*k+mm*k+r] = (n*lambda*pow(pi[r],gamma) - S[j*m*k+mm*k+r]) - / Gram2[j*p*k+j*k+r]; + S[ai(j,mm,r,p,m,k)] = -rho[ai(mm,mm,r,m,m,k)] * ps2[ai(j,mm,r,p,m,k)] + dotPhiGram2; + if (fabs(S[ai(j,mm,r,p,m,k)]) <= n*lambda*pow(pi[r],gamma)) + phi[ai(j,mm,r,p,m,k)] = 0; + else if (S[ai(j,mm,r,p,m,k)] > n*lambda*pow(pi[r],gamma)) + phi[ai(j,mm,r,p,m,k)] = (n*lambda*pow(pi[r],gamma) - S[ai(j,mm,r,p,m,k)]) + / Gram2[ai(j,j,r,p,p,k)]; else - phi[j*m*k+mm*k+r] = -(n*lambda*pow(pi[r],gamma) + S[j*m*k+mm*k+r]) - / Gram2[j*p*k+j*k+r]; + phi[ai(j,mm,r,p,m,k)] = -(n*lambda*pow(pi[r],gamma) + S[ai(j,mm,r,p,m,k)]) + / Gram2[ai(j,j,r,p,p,k)]; } } } - + ///////////// // Etape E // ///////////// - + int signum; - Real sumLogLLF2 = 0.0; - for (mwSize i=0; i<n; i++) + double sumLogLLF2 = 0.0; + for (int i=0; i<n; i++) { - Real sumLLF1 = 0.0; - Real sumGamI = 0.0; - Real minDotProduct = INFINITY; - - for (mwSize r=0; r<k; r++) + double sumLLF1 = 0.0; + double sumGamI = 0.0; + double minDotProduct = INFINITY; + + for (int r=0; r<k; r++) { //Compute //Gam(i,r) = Pi(r) * det(Rho(:,:,r)) * exp( -1/2 * (Y(i,:)*Rho(:,:,r) - X(i,:)... @@ -274,57 +274,57 @@ void EMGLLF( //split in several sub-steps //compute Y(i,:)*rho(:,:,r) - for (mwSize u=0; u<m; u++) + for (int u=0; u<m; u++) { YiRhoR[u] = 0.0; - for (mwSize v=0; v<m; v++) - YiRhoR[u] += Y[i*m+v] * rho[v*m*k+u*k+r]; + for (int v=0; v<m; v++) + YiRhoR[u] += Y[imi(i,v,n,m)] * rho[ai(v,u,r,m,m,k)]; } - + //compute X(i,:)*phi(:,:,r) - for (mwSize u=0; u<m; u++) + for (int u=0; u<m; u++) { XiPhiR[u] = 0.0; - for (mwSize v=0; v<p; v++) - XiPhiR[u] += X[i*p+v] * phi[v*m*k+u*k+r]; + for (int v=0; v<p; v++) + XiPhiR[u] += X[mi(i,v,n,p)] * phi[ai(v,u,r,p,m,k)]; } - + // compute dotProduct < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) > dotProducts[r] = 0.0; - for (mwSize u=0; u<m; u++) + for (int u=0; u<m; u++) dotProducts[r] += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]); if (dotProducts[r] < minDotProduct) minDotProduct = dotProducts[r]; } - Real shift = 0.5*minDotProduct; - for (mwSize r=0; r<k; r++) + double shift = 0.5*minDotProduct; + for (int r=0; r<k; r++) { //compute det(rho(:,:,r)) [TODO: avoid re-computations] - for (mwSize u=0; u<m; u++) + for (int u=0; u<m; u++) { - for (mwSize v=0; v<m; v++) - matrix->data[u*m+v] = rho[u*m*k+v*k+r]; + for (int v=0; v<m; v++) + matrix->data[u*m+v] = rho[ai(u,v,r,m,m,k)]; } gsl_linalg_LU_decomp(matrix, permutation, &signum); - Real detRhoR = gsl_linalg_LU_det(matrix, signum); - - Gam[i*k+r] = pi[r] * detRhoR * exp(-0.5*dotProducts[r] + shift); - sumLLF1 += Gam[i*k+r] / pow(2*M_PI,m/2.0); - sumGamI += Gam[i*k+r]; + double detRhoR = gsl_linalg_LU_det(matrix, signum); + + Gam[mi(i,r,n,k)] = pi[r] * detRhoR * exp(-0.5*dotProducts[r] + shift); + sumLLF1 += Gam[mi(i,r,n,k)] / pow(2*M_PI,m/2.0); + sumGamI += Gam[mi(i,r,n,k)]; } sumLogLLF2 += log(sumLLF1); - for (mwSize r=0; r<k; r++) + for (int r=0; r<k; r++) { //gam(i,r)=Gam(i,r)/sum(Gam(i,:)); - gam[i*k+r] = sumGamI > EPS - ? Gam[i*k+r] / sumGamI + gam[mi(i,r,n,k)] = sumGamI > EPS + ? Gam[mi(i,r,n,k)] / sumGamI : 0.0; } } //sum(pen(ite,:)) - Real sumPen = 0.0; - for (mwSize r=0; r<k; r++) + double sumPen = 0.0; + for (int r=0; r<k; r++) sumPen += pow(pi[r],gamma) * b[r]; //LLF(ite)=-1/n*sum(log(LLF2(ite,:)))+lambda*sum(pen(ite,:)); LLF[ite] = -invN * sumLogLLF2 + lambda * sumPen; @@ -334,42 +334,42 @@ void EMGLLF( dist = (LLF[ite] - LLF[ite-1]) / (1.0 + fabs(LLF[ite])); //Dist1=max(max((abs(phi-Phi))./(1+abs(phi)))); - Real Dist1 = 0.0; - for (mwSize u=0; u<p; u++) + double Dist1 = 0.0; + for (int u=0; u<p; u++) { - for (mwSize v=0; v<m; v++) + for (int v=0; v<m; v++) { - for (mwSize w=0; w<k; w++) + for (int w=0; w<k; w++) { - Real tmpDist = fabs(phi[u*m*k+v*k+w]-Phi[u*m*k+v*k+w]) - / (1.0+fabs(phi[u*m*k+v*k+w])); + double tmpDist = fabs(phi[ai(u,v,w,p,m,k)]-Phi[ai(u,v,w,p,m,k)]) + / (1.0+fabs(phi[ai(u,v,w,p,m,k)])); if (tmpDist > Dist1) Dist1 = tmpDist; } } } //Dist2=max(max((abs(rho-Rho))./(1+abs(rho)))); - Real Dist2 = 0.0; - for (mwSize u=0; u<m; u++) + double Dist2 = 0.0; + for (int u=0; u<m; u++) { - for (mwSize v=0; v<m; v++) + for (int v=0; v<m; v++) { - for (mwSize w=0; w<k; w++) + for (int w=0; w<k; w++) { - Real tmpDist = fabs(rho[u*m*k+v*k+w]-Rho[u*m*k+v*k+w]) - / (1.0+fabs(rho[u*m*k+v*k+w])); + double tmpDist = fabs(rho[ai(u,v,w,m,m,k)]-Rho[ai(u,v,w,m,m,k)]) + / (1.0+fabs(rho[ai(u,v,w,m,m,k)])); if (tmpDist > Dist2) Dist2 = tmpDist; } } } //Dist3=max(max((abs(pi-Pi))./(1+abs(Pi)))); - Real Dist3 = 0.0; - for (mwSize u=0; u<n; u++) + double Dist3 = 0.0; + for (int u=0; u<n; u++) { - for (mwSize v=0; v<k; v++) + for (int v=0; v<k; v++) { - Real tmpDist = fabs(pi[v]-Pi[v]) / (1.0+fabs(pi[v])); + double tmpDist = fabs(pi[v]-Pi[v]) / (1.0+fabs(pi[v])); if (tmpDist > Dist3) Dist3 = tmpDist; } diff --git a/src/sources/utils/io.c b/src/sources/utils/io.c deleted file mode 100644 index 36a1d8e..0000000 --- a/src/sources/utils/io.c +++ /dev/null @@ -1,168 +0,0 @@ -#include "ioutils.h" -#include <string.h> -#include <stdio.h> - -// Check if array == refArray -void compareArray(const char* ID, const void* array, const void* refArray, mwSize size, int isInteger) -{ - Real EPS = 1e-5; //precision - printf("Checking %s\n",ID); - Real maxError = 0.0; - for (mwSize i=0; i<size; i++) - { - Real error = isInteger - ? fabs(((Int*)array)[i] - ((Int*)refArray)[i]) - : fabs(((Real*)array)[i] - ((Real*)refArray)[i]); - if (error >= maxError) - maxError = error; - } - if (maxError >= EPS) - printf(" Inaccuracy: max(abs(error)) = %g >= %g\n",maxError,EPS); - else - printf(" OK\n"); -} - -// Next function is a generalization of : -//~ Real* brToMatlabArray(Real* brArray, int dimX, int dimY, int dimZ, int dimW) -//~ { - //~ Real* matlabArray = (Real*)malloc(dimX*dimY*dimZ*dimW*sizeof(Real)); - //~ for (int u=0; u<dimX*dimY*dimZ*dimW; u++) - //~ { - //~ int xIndex = u / (dimY*dimZ*dimW); - //~ int yIndex = (u % (dimY*dimZ*dimW)) / (dimZ*dimW); - //~ int zIndex = (u % (dimZ*dimW)) / dimW; - //~ int wIndex = u % dimW; - //~ matlabArray[xIndex+yIndex*dimX+zIndex*dimX*dimY+wIndex*dimX*dimY*dimZ] = brArray[u]; - //~ } - //~ return matlabArray; -//~ } - -// Auxiliary to convert from ours ("by-rows") encoding to MATLAB -void* brToMatlabArray(const void* brArray, const mwSize* dimensions, int nbDims, int isInteger) -{ - mwSize totalDim = 1; - for (int i=0; i<nbDims; i++) - totalDim *= dimensions[i]; - size_t elementSize = isInteger - ? sizeof(Int) - : sizeof(Real); - - void* matlabArray = malloc(totalDim*elementSize); - for (mwSize u=0; u<totalDim; u++) - { - mwSize prodDimLeft = totalDim; - mwSize prodDimRight = totalDim / dimensions[0]; - mwSize prodDimInIndex = 1; - mwSize index = 0; - for (int v=0; v<nbDims; v++) - { - index += ((u % prodDimLeft) / prodDimRight) * prodDimInIndex; - prodDimInIndex *= dimensions[v]; - prodDimLeft /= dimensions[v]; - if (v+1 < nbDims) - prodDimRight /= dimensions[v+1]; - } - if (isInteger) - ((Int*)matlabArray)[index] = ((Int*)brArray)[u]; - else - ((Real*)matlabArray)[index] = ((Real*)brArray)[u]; - } - return matlabArray; -} - -// Next function is a generalization of : -//~ Real* matlabToBrArray(Real* matlabArray, int dimX, int dimY, int dimZ, int dimU) -//~ { - //~ Real* brArray = (Real*)malloc(dimX*dimY*dimZ*dimU*sizeof(Real)); - //~ for (int u=0; u<dimX*dimY*dimZ*dimU; u++) - //~ { - //~ int xIndex = u % dimX; - //~ int yIndex = (u % (dimX*dimY)) / dimX; - //~ int zIndex = (u % (dimX*dimY*dimZ)) / (dimX*dimY); - //~ int uIndex = u / (dimX*dimY*dimZ); - //~ brArray[xIndex*dimY*dimZ*dimU+yIndex*dimZ*dimU+zIndex*dimU+uIndex] = matlabArray[u]; - //~ } - //~ return brArray; -//~ } - -// Auxiliary to convert from MATLAB encoding to ours ("by-rows") -void* matlabToBrArray(const void* matlabArray, const mwSize* dimensions, int nbDims, int isInteger) -{ - mwSize totalDim = 1; - for (int i=0; i<nbDims; i++) - totalDim *= dimensions[i]; - size_t elementSize = isInteger - ? sizeof(Int) - : sizeof(Real); - - void* brArray = malloc(totalDim*elementSize); - for (mwSize u=0; u<totalDim; u++) - { - mwSize prodDimLeft = dimensions[0]; - mwSize prodDimRight = 1; - mwSize prodDimInIndex = totalDim / dimensions[0]; - mwSize index = 0; - for (int v=0; v<nbDims; v++) - { - index += ((u % prodDimLeft) / prodDimRight) * prodDimInIndex; - if (v+1 < nbDims) - { - prodDimInIndex /= dimensions[v+1]; - prodDimLeft *= dimensions[v+1]; - } - prodDimRight *= dimensions[v]; - } - if (isInteger) - ((Int*)brArray)[index] = ((Int*)matlabArray)[u]; - else - ((Real*)brArray)[index] = ((Real*)matlabArray)[u]; - } - return brArray; -} - -// Read array by columns (as in MATLAB) and return by-rows encoding -void* readArray(const char* fileName, const mwSize* dimensions, int nbDims, int isInteger) -{ - // need to prepend '../data/' (not really nice code...) - char* fullFileName = (char*)calloc(8+strlen(fileName)+1,sizeof(char)); - strcat(fullFileName, "../data/"); - strcat(fullFileName, fileName); - FILE* file = fopen(fullFileName, "r"); - free(fullFileName); - - mwSize totalDim = 1; - for (int i=0; i<nbDims; i++) - totalDim *= dimensions[i]; - size_t elementSize = isInteger - ? sizeof(Int) - : sizeof(Real); - - // read all values, and convert them to by-rows matrices format - void* matlabArray = malloc(totalDim*elementSize); - char curChar = ' '; - char bufferNum[64]; - for (mwSize u=0; u<totalDim; u++) - { - // position to next non-separator character - while (!feof(file) && (curChar==' ' || curChar=='\n' || curChar=='\t' || curChar==',')) - curChar = fgetc(file); - // read number (as a string) - int bufferIndex = 0; - while (!feof(file) && curChar!=' ' && curChar!='\n' && curChar!='\t' && curChar!=',') - { - bufferNum[bufferIndex++] = curChar; - curChar = fgetc(file); - } - bufferNum[bufferIndex] = 0; - // transform string into Real, and store it at appropriate location - if (isInteger) - ((Int*)matlabArray)[u] = atoi(bufferNum); - else - ((Real*)matlabArray)[u] = atof(bufferNum); - } - fclose(file); - - void* brArray = matlabToBrArray(matlabArray, dimensions, nbDims, isInteger); - free(matlabArray); - return brArray; -} diff --git a/src/sources/utils/io.h b/src/sources/utils/io.h index 31f464a..9aa5899 100644 --- a/src/sources/utils/io.h +++ b/src/sources/utils/io.h @@ -1,75 +1,26 @@ #ifndef select_ioutils_H #define select_ioutils_H -#include <stdlib.h> -#include <math.h> -#include <stdint.h> -#include <uchar.h> //for type wchar16_t - -// Include header for mwSize type -#ifdef Octave -#include <mex.h> -#else -#include <tmwtypes.h> -#endif - -// CHUNK_SIZE = number of lambda values to be treated sequentially by a single core -#define CHUNK_SIZE 1 - -// integer type chosen in MATLAB (to be tuned: 32 bits should be enough) -typedef int64_t Int; - -// real number type chosen in MATLAB (default: double) -typedef double Real; - -#ifndef M_PI -#define M_PI 3.141592653589793 -#endif - // Fill an array with zeros #define zeroArray(array, size)\ {\ - for (Int u=0; u<size; u++)\ + for (int u=0; u<size; u++)\ array[u] = 0;\ } // Copy an 1D array #define copyArray(array, copy, size)\ {\ - for (Int u=0; u<size; u++)\ + for (int u=0; u<size; u++)\ copy[u] = array[u];\ } -// Check if array == refArray -void compareArray(const char* ID, const void* array, const void* refArray, mwSize size, int isInteger); - -#define compareArray_int(ID, array, refArray, size)\ - compareArray(ID, array, refArray, size, 1) -#define compareArray_real(ID, array, refArray, size)\ - compareArray(ID, array, refArray, size, 0) - -// Auxiliary to convert from ours ("by-rows") encoding to MATLAB -void* brToMatlabArray(const void* brArray, const mwSize* dimensions, int nbDims, int isInteger); - -#define brToMatlabArray_int(brArray, dimensions, nbDims)\ - (Int*)brToMatlabArray(brArray, dimensions, nbDims, 1) -#define brToMatlabArray_real(brArray, dimensions, nbDims)\ - (Real*)brToMatlabArray(brArray, dimensions, nbDims, 0) - -// Auxiliary to convert from MATLAB encoding to ours ("by-rows") -void* matlabToBrArray(const void* matlabArray, const mwSize* dimensions, int nbDims, int isInteger); - -#define matlabToBrArray_int(matlabArray, dimensions, nbDims)\ - (Int*)matlabToBrArray(matlabArray, dimensions, nbDims, 1) -#define matlabToBrArray_real(matlabArray, dimensions, nbDims)\ - (Real*)matlabToBrArray(matlabArray, dimensions, nbDims, 0) - -// Read array by columns (as in MATLAB) and return by-rows encoding -void* readArray(const char* fileName, const mwSize* dimensions, int nbDims, int isInteger); +// Matrix Index ; TODO? ncol unused +#define mi(i,j,nrow,ncol)\ + j*nrow + i -#define readArray_int(fileName, dimensions, nbDims)\ - (Int*)readArray(fileName, dimensions, nbDims, 1) -#define readArray_real(fileName, dimensions, nbDims)\ - (Real*)readArray(fileName, dimensions, nbDims, 0) +// Array Index ; TODO? d3 unused +#define ai(i,j,k,d1,d2,d3)\ + k*d1*d2 + j*d1 + i #endif diff --git a/src/sources/utils/omp_num_threads.h b/src/sources/utils/omp_num_threads.h deleted file mode 100644 index f4fcc24..0000000 --- a/src/sources/utils/omp_num_threads.h +++ /dev/null @@ -1,7 +0,0 @@ -#ifndef select_omp_num_threads_H -#define select_omp_num_threads_H - -//Par exemple... -#define OMP_NUM_THREADS 8 - -#endif diff --git a/src/sources/utils/tune_parallelisms.h b/src/sources/utils/tune_parallelisms.h new file mode 100644 index 0000000..33a5171 --- /dev/null +++ b/src/sources/utils/tune_parallelisms.h @@ -0,0 +1,10 @@ +#ifndef tune_parallelism_H +#define tune_parallelism_H + +// Number of OpenMP threads +#define OMP_NUM_THREADS 8 + +// CHUNK_SIZE = number of lambda values to be treated sequentially by a single core +#define CHUNK_SIZE 1 + +#endif -- 2.44.0