From: Benjamin Auder Date: Fri, 2 Dec 2016 01:41:42 +0000 (+0100) Subject: finished a.EMGLLF.c [to test] ; fix indices in EMGLLF.c [by cols] X-Git-Url: https://git.auder.net/js/css/%7B%7B?a=commitdiff_plain;h=4cab944a7bf905e811d740694c39f16b2b23d7e1;p=valse.git finished a.EMGLLF.c [to test] ; fix indices in EMGLLF.c [by cols] --- 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 #include #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 -// 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 - Real dotProduct = 0.0; - for (mwSize u=0; u 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 dotProducts[r] = 0.0; - for (mwSize u=0; udata[u*m+v] = rho[u*m*k+v*k+r]; + for (int v=0; vdata[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 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 Dist1) Dist1 = tmpDist; } } } //Dist2=max(max((abs(rho-Rho))./(1+abs(rho)))); - Real Dist2 = 0.0; - for (mwSize u=0; u Dist2) Dist2 = tmpDist; } } } //Dist3=max(max((abs(pi-Pi))./(1+abs(Pi)))); - Real Dist3 = 0.0; - for (mwSize u=0; u 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 -#include - -// 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= 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 -#include -#include -#include //for type wchar16_t - -// Include header for mwSize type -#ifdef Octave -#include -#else -#include -#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