X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=src%2Fsources%2FEMGrank.c;h=3a9bf94b23a05c443d7cb6f32dd835d19667274f;hb=53fa233d8fbeaf4d51a4874ba69d8472d01d04ba;hp=2e85ce6fd019f9c836cc26d82210b1fc7d3a1556;hpb=552b00e200e8a990c1247989d29e98d4ae8679f3;p=valse.git

diff --git a/src/sources/EMGrank.c b/src/sources/EMGrank.c
index 2e85ce6..3a9bf94 100644
--- a/src/sources/EMGrank.c
+++ b/src/sources/EMGrank.c
@@ -1,14 +1,15 @@
-#include "EMGrank.h"
+#include <stdlib.h>
 #include <gsl/gsl_linalg.h>
+#include "utils.h"
 
 // Compute pseudo-inverse of a square matrix
-static double* pinv(const double* matrix, int dim)
+static Real* pinv(const Real* 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"
+	Real 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));
+	Real* inverse = (Real*)malloc(dim*dim*sizeof(Real));
 	for (int i=0; i<dim; i++)
 	{
 		for (int ii=0; ii<dim; ii++)
 		{
-			double dotProduct = 0.0;
+			Real dotProduct = 0.0;
 			for (int j=0; j<dim; j++)
 				dotProduct += V->data[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 Real* Pi, // parametre de proportion
+	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 double* X, // régresseurs
-	const double* Y, // réponse
-	double tau, // seuil pour accepter la convergence
+	const Real* X, // régresseurs
+	const Real* Y, // réponse
+	Real 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
+	Real* phi, // parametre de moyenne renormalisé, calculé par l'EM
+	Real* 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));
+	Real* Phi = (Real*)calloc(p*m*k,sizeof(Real));
+	Real* hatBetaR = (Real*)malloc(p*m*sizeof(Real));
 	int signum;
-	double invN = 1.0/n;
+	Real invN = 1.0/n;
 	int deltaPhiBufferSize = 20;
-	double* deltaPhi = (double*)malloc(deltaPhiBufferSize*sizeof(double));
+	Real* deltaPhi = (Real*)malloc(deltaPhiBufferSize*sizeof(Real));
 	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));
+	Real sumDeltaPhi = 0.0;
+	Real* YiRhoR = (Real*)malloc(m*sizeof(Real));
+	Real* XiPhiR = (Real*)malloc(m*sizeof(Real));
+	Real* Xr = (Real*)malloc(n*p*sizeof(Real));
+	Real* Yr = (Real*)malloc(n*m*sizeof(Real));
+	Real* tXrXr = (Real*)malloc(p*p*sizeof(Real));
+	Real* tXrYr = (Real*)malloc(p*m*sizeof(Real));
 	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; i<p*m*k; i++)
-		phi[i] = 0.0;
+	zeroArray(phi, p*m*k);
 
 	while (ite<mini || (ite<maxi && sumDeltaPhi>tau))
 	{
@@ -115,7 +115,7 @@ void EMGrank(
 			{
 				for (int jj=0; jj<p; jj++)
 				{
-					double dotProduct = 0.0;
+					Real dotProduct = 0.0;
 					for (int u=0; u<cardClustR; u++)
 						dotProduct += Xr[mi(u,j,n,p)] * Xr[mi(u,jj,n,p)];
 					tXrXr[mi(j,jj,p,p)] = dotProduct;
@@ -123,16 +123,16 @@ void EMGrank(
 			}
 
 			//Get pseudo inverse = (t(Xr)*Xr)^{-1}
-			double* invTXrXr = pinv(tXrXr, p);
+			Real* invTXrXr = pinv(tXrXr, p);
 			
 			// Compute tXrYr = t(Xr) * Yr
 			for (int j=0; j<p; j++)
 			{
 				for (int jj=0; jj<m; jj++)
 				{
-					double dotProduct = 0.0;
+					Real dotProduct = 0.0;
 					for (int u=0; u<cardClustR; u++)
-						dotProduct += Xr[mi(u,j,n,p)] * Yr[mi(u,j,n,m)];
+						dotProduct += Xr[mi(u,j,n,p)] * Yr[mi(u,jj,n,m)];
 					tXrYr[mi(j,jj,p,m)] = dotProduct;
 				}
 			}
@@ -142,7 +142,7 @@ void EMGrank(
 			{
 				for (int jj=0; jj<m; jj++)
 				{
-					double dotProduct = 0.0;
+					Real dotProduct = 0.0;
 					for (int u=0; u<p; u++)
 						dotProduct += invTXrXr[mi(j,u,p,p)] * tXrYr[mi(u,jj,p,m)];
 					matrixM->data[j*m+jj] = dotProduct;
@@ -159,12 +159,12 @@ void EMGrank(
 				S->data[j] = 0.0;
 			
 			//[intermediate step] Compute hatBetaR = U * S * t(V)
-			double* U = matrixM->data;
+			double* U = matrixM->data; //GSL require double precision
 			for (int j=0; j<p; j++)
 			{
 				for (int jj=0; jj<m; jj++)
 				{
-					double dotProduct = 0.0;
+					Real dotProduct = 0.0;
 					for (int u=0; u<m; u++)
 						dotProduct += U[j*m+u] * S->data[u] * V->data[jj*m+u];
 					hatBetaR[mi(j,jj,p,m)] = dotProduct;
@@ -176,7 +176,7 @@ void EMGrank(
 			{
 				for (int jj=0; jj<m; jj++)
 				{
-					double dotProduct=0.0;
+					Real dotProduct=0.0;
 					for (int u=0; u<m; u++)
 						dotProduct += hatBetaR[mi(j,u,p,m)] * Rho[ai(u,jj,r,m,m,k)];
 					phi[ai(j,jj,r,p,m,k)] = dotProduct;
@@ -188,11 +188,11 @@ void EMGrank(
 		// Etape E //
 		/////////////
 		
-		double sumLogLLF2 = 0.0;
+		Real sumLogLLF2 = 0.0;
 		for (int i=0; i<n; i++)
 		{
-			double sumLLF1 = 0.0;
-			double maxLogGamIR = -INFINITY;
+			Real sumLLF1 = 0.0;
+			Real maxLogGamIR = -INFINITY;
 			for (int r=0; r<k; r++)
 			{
 				//Compute
@@ -207,7 +207,7 @@ void EMGrank(
 						matrixE->data[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);
+				Real detRhoR = gsl_linalg_LU_det(matrixE, signum);
 
 				//compute Y(i,:)*Rho(:,:,r)
 				for (int j=0; j<m; j++)
@@ -226,10 +226,10 @@ void EMGrank(
 				}
 
 				//compute dotProduct < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) >
-				double dotProduct = 0.0;
+				Real dotProduct = 0.0;
 				for (int u=0; u<m; u++)
 					dotProduct += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
-				double logGamIR = log(Pi[r]) + log(detRhoR) - 0.5*dotProduct;
+				Real logGamIR = log(Pi[r]) + log(detRhoR) - 0.5*dotProduct;
 
 				//Z(i) = index of max (gam(i,:))
 				if (logGamIR > maxLogGamIR)
@@ -247,14 +247,14 @@ void EMGrank(
 		*LLF = -invN * sumLogLLF2;
 
 		//newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi))));
-		double newDeltaPhi = 0.0;
+		Real newDeltaPhi = 0.0;
 		for (int j=0; j<p; j++)
 		{
 			for (int jj=0; jj<m; jj++)
 			{
 				for (int r=0; r<k; r++)
 				{
-					double tmpDist = fabs(phi[ai(j,jj,r,p,m,k)]-Phi[ai(j,jj,r,p,m,k)])
+					Real tmpDist = fabs(phi[ai(j,jj,r,p,m,k)]-Phi[ai(j,jj,r,p,m,k)])
 						/ (1.0+fabs(phi[ai(j,jj,r,p,m,k)]));
 					if (tmpDist > newDeltaPhi)
 						newDeltaPhi = tmpDist;
@@ -281,7 +281,7 @@ void EMGrank(
 			for (int jj=0; jj<p; jj++)
 			{
 				for (int r=0; r<k; r++)
-					Phi[ai(j,jj,r,p,m,k)] = phi[ai(j,jj,r,p,m,k)];
+					Phi[ai(jj,j,r,p,m,k)] = phi[ai(jj,j,r,p,m,k)];
 			}
 		}
 		ite++;