X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=src%2Fsources%2FEMGLLF.c;h=087116be27573bc679b8b41c0d264ae0ece95700;hb=9ff729fb6afff5bed327fa8619138fd3b6f6f13b;hp=0c39d6e05513da54c88c4a24f2918b529c7b7654;hpb=4cab944a7bf905e811d740694c39f16b2b23d7e1;p=valse.git diff --git a/src/sources/EMGLLF.c b/src/sources/EMGLLF.c index 0c39d6e..087116b 100644 --- a/src/sources/EMGLLF.c +++ b/src/sources/EMGLLF.c @@ -1,31 +1,32 @@ -#include "EMGLLF.h" +#include "utils.h" +#include #include // TODO: don't recompute indexes every time...... -void EMGLLF( +void EMGLLF_core( // IN parameters - const double* phiInit, // parametre initial de moyenne renormalisé - const double* rhoInit, // parametre initial de variance renormalisé - 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 + 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, // 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 // OUT parameters (all pointers, to be modified) - 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, + 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ée à cet échantillon, pour les valeurs estimées des paramètres + Real* S, // additional size parameters - 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 + 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); @@ -36,32 +37,32 @@ void EMGLLF( //Other local variables //NOTE: variables order is always [maxi],n,p,m,k - double* gam = (double*)malloc(n*k*sizeof(double)); + Real* gam = (Real*)malloc(n*k*sizeof(Real)); copyArray(gamInit, gam, n*k); - 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)); + 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)); gsl_matrix* matrix = gsl_matrix_alloc(m, m); gsl_permutation* permutation = gsl_permutation_alloc(m); - double* YiRhoR = (double*)malloc(m*sizeof(double)); - double* XiPhiR = (double*)malloc(m*sizeof(double)); - double dist = 0.; - double dist2 = 0.; + Real* YiRhoR = (Real*)malloc(m*sizeof(Real)); + Real* XiPhiR = (Real*)malloc(m*sizeof(Real)); + Real dist = 0.; + Real dist2 = 0.; int ite = 0; - double EPS = 1e-15; - double* dotProducts = (double*)malloc(k*sizeof(double)); + Real EPS = 1e-15; + Real* dotProducts = (Real*)malloc(k*sizeof(Real)); while (ite < mini || (ite < maxi && (dist >= tau || dist2 >= sqrt(tau)))) { @@ -89,7 +90,7 @@ void EMGLLF( //ps2(:,mm,r)=transpose(X2(:,:,r))*Y2(:,mm,r); for (int u=0; u - double dotProduct = 0.0; + Real dotProduct = 0.0; for (int u=0; u 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)]) + 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[ai(j,mm,r,p,m,k)] = -(n*lambda*pow(pi[r],gamma) + S[ai(j,mm,r,p,m,k)]) + 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)]; } } @@ -259,18 +262,18 @@ void EMGLLF( ///////////// int signum; - double sumLogLLF2 = 0.0; + Real sumLogLLF2 = 0.0; for (int i=0; i + //compute dotProduct + // < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) > dotProducts[r] = 0.0; for (int u=0; udata[u*m+v] = rho[ai(u,v,r,m,m,k)]; } gsl_linalg_LU_decomp(matrix, permutation, &signum); - double detRhoR = gsl_linalg_LU_det(matrix, signum); + Real 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); @@ -323,7 +327,7 @@ void EMGLLF( } //sum(pen(ite,:)) - double sumPen = 0.0; + Real sumPen = 0.0; for (int r=0; r Dist1) Dist1 = tmpDist; @@ -349,14 +353,14 @@ void EMGLLF( } } //Dist2=max(max((abs(rho-Rho))./(1+abs(rho)))); - double Dist2 = 0.0; + Real Dist2 = 0.0; for (int u=0; u Dist2) Dist2 = tmpDist; @@ -364,12 +368,12 @@ void EMGLLF( } } //Dist3=max(max((abs(pi-Pi))./(1+abs(Pi)))); - double Dist3 = 0.0; + Real Dist3 = 0.0; for (int u=0; u Dist3) Dist3 = tmpDist; }