X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=src%2Fsources%2FEMGrank.c;h=3d5a9c760b807f978033b68ff3b661e22e861662;hb=afa07d41c7592ac0ccd55d7af23c3bfef213291e;hp=2e85ce6fd019f9c836cc26d82210b1fc7d3a1556;hpb=552b00e200e8a990c1247989d29e98d4ae8679f3;p=valse.git diff --git a/src/sources/EMGrank.c b/src/sources/EMGrank.c index 2e85ce6..3d5a9c7 100644 --- a/src/sources/EMGrank.c +++ b/src/sources/EMGrank.c @@ -1,14 +1,15 @@ -#include "EMGrank.h" +#include #include +#include "utils.h" // Compute pseudo-inverse of a square matrix -static double* pinv(const double* matrix, int dim) +static float* pinv(const float* matrix, int dim) { gsl_matrix* U = gsl_matrix_alloc(dim,dim); gsl_matrix* V = gsl_matrix_alloc(dim,dim); gsl_vector* S = gsl_vector_alloc(dim); gsl_vector* work = gsl_vector_alloc(dim); - double EPS = 1e-10; //threshold for singular value "== 0" + float EPS = 1e-10; //threshold for singular value "== 0" //copy matrix into U copyArray(matrix, U->data, dim*dim); @@ -18,12 +19,12 @@ static double* pinv(const double* matrix, int dim) gsl_vector_free(work); // Obtain pseudo-inverse by V*S^{-1}*t(U) - double* inverse = (double*)malloc(dim*dim*sizeof(double)); + float* inverse = (float*)malloc(dim*dim*sizeof(float)); for (int i=0; idata[i*dim+j] * (S->data[j] > EPS ? 1.0/S->data[j] : 0.0) * U->data[ii*dim+j]; inverse[i*dim+ii] = dotProduct; @@ -37,19 +38,19 @@ static double* pinv(const double* matrix, int dim) } // TODO: comment EMGrank purpose -void EMGrank( +void EMGrank_core( // IN parameters - const double* Pi, // parametre de proportion - const double* Rho, // parametre initial de variance renormalisé + const float* Pi, // parametre de proportion + const float* 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 double* X, // régresseurs - const double* Y, // réponse - double tau, // seuil pour accepter la convergence + const float* X, // régresseurs + const float* Y, // réponse + float tau, // seuil pour accepter la convergence const int* rank, // vecteur des rangs possibles // OUT parameters - double* phi, // parametre de moyenne renormalisé, calculé par l'EM - double* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres + float* phi, // parametre de moyenne renormalisé, calculé par l'EM + float* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres // additional size parameters int n, // taille de l'echantillon int p, // nombre de covariables @@ -57,20 +58,20 @@ void EMGrank( int k) // nombre de composantes { // Allocations, initializations - double* Phi = (double*)calloc(p*m*k,sizeof(double)); - double* hatBetaR = (double*)malloc(p*m*sizeof(double)); + float* Phi = (float*)calloc(p*m*k,sizeof(float)); + float* hatBetaR = (float*)malloc(p*m*sizeof(float)); int signum; - double invN = 1.0/n; + float invN = 1.0/n; int deltaPhiBufferSize = 20; - double* deltaPhi = (double*)malloc(deltaPhiBufferSize*sizeof(double)); + float* deltaPhi = (float*)malloc(deltaPhiBufferSize*sizeof(float)); int ite = 0; - double sumDeltaPhi = 0.0; - double* YiRhoR = (double*)malloc(m*sizeof(double)); - double* XiPhiR = (double*)malloc(m*sizeof(double)); - double* Xr = (double*)malloc(n*p*sizeof(double)); - double* Yr = (double*)malloc(n*m*sizeof(double)); - double* tXrXr = (double*)malloc(p*p*sizeof(double)); - double* tXrYr = (double*)malloc(p*m*sizeof(double)); + float sumDeltaPhi = 0.0; + float* YiRhoR = (float*)malloc(m*sizeof(float)); + float* XiPhiR = (float*)malloc(m*sizeof(float)); + float* Xr = (float*)malloc(n*p*sizeof(float)); + float* Yr = (float*)malloc(n*m*sizeof(float)); + float* tXrXr = (float*)malloc(p*p*sizeof(float)); + float* tXrYr = (float*)malloc(p*m*sizeof(float)); gsl_matrix* matrixM = gsl_matrix_alloc(p, m); gsl_matrix* matrixE = gsl_matrix_alloc(m, m); gsl_permutation* permutation = gsl_permutation_alloc(m); @@ -82,8 +83,7 @@ void EMGrank( int* Z = (int*)calloc(n, sizeof(int)); //Initialize phi to zero, because some M loops might exit before phi affectation - for (int i=0; itau)) { @@ -115,7 +115,7 @@ void EMGrank( { for (int jj=0; jjdata[j*m+jj] = dotProduct; @@ -164,7 +164,7 @@ void EMGrank( { for (int jj=0; jjdata[u] * V->data[jj*m+u]; hatBetaR[mi(j,jj,p,m)] = dotProduct; @@ -176,7 +176,7 @@ void EMGrank( { for (int jj=0; jjdata[j*m+jj] = Rho[ai(j,jj,r,m,m,k)]; } gsl_linalg_LU_decomp(matrixE, permutation, &signum); - double detRhoR = gsl_linalg_LU_det(matrixE, signum); + float detRhoR = gsl_linalg_LU_det(matrixE, signum); //compute Y(i,:)*Rho(:,:,r) for (int j=0; j - double dotProduct = 0.0; + float dotProduct = 0.0; for (int u=0; u maxLogGamIR) @@ -247,14 +247,14 @@ void EMGrank( *LLF = -invN * sumLogLLF2; //newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi)))); - double newDeltaPhi = 0.0; + float newDeltaPhi = 0.0; for (int j=0; j newDeltaPhi) newDeltaPhi = tmpDist;