[valse.git] / src / sources / EMGrank.c

#include "EMGrank.h"
#include <gsl/gsl_linalg.h>

// Compute pseudo-inverse of a square matrix
static double* pinv(const double* 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"
	
	//copy matrix into U
	copyArray(matrix, U->data, dim*dim);

	//U,S,V = SVD of matrix
	gsl_linalg_SV_decomp(U, V, S, work);
	gsl_vector_free(work);

	// Obtain pseudo-inverse by V*S^{-1}*t(U)
	double* inverse = (double*)malloc(dim*dim*sizeof(double));
	for (int i=0; i<dim; i++)
	{
		for (int ii=0; ii<dim; ii++)
		{
			double 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;
		}
	}
	
	gsl_matrix_free(U);
	gsl_matrix_free(V);
	gsl_vector_free(S);
	return inverse;
}

// TODO: comment EMGrank purpose
void EMGrank(
	// IN parameters
	const double* Pi, // parametre de proportion
	const double* 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 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
	// additional size parameters
	int n, // taille de l'echantillon
	int p, // nombre de covariables
	int m, // taille de Y (multivarié)
	int k) // nombre de composantes
{
	// Allocations, initializations
	double* Phi = (double*)calloc(p*m*k,sizeof(double));
	double* hatBetaR = (double*)malloc(p*m*sizeof(double));
	int signum;
	double invN = 1.0/n;
	int deltaPhiBufferSize = 20;
	double* deltaPhi = (double*)malloc(deltaPhiBufferSize*sizeof(double));
	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));
	gsl_matrix* matrixM = gsl_matrix_alloc(p, m);
	gsl_matrix* matrixE = gsl_matrix_alloc(m, m);
	gsl_permutation* permutation = gsl_permutation_alloc(m);
	gsl_matrix* V = gsl_matrix_alloc(m,m);
	gsl_vector* S = gsl_vector_alloc(m);
	gsl_vector* work = gsl_vector_alloc(m);

	//Initialize class memberships (all elements in class 0; TODO: randomize ?)
	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;

	while (ite<mini || (ite<maxi && sumDeltaPhi>tau))
	{
		/////////////
		// Etape M //
		/////////////
		
		//M step: Mise à jour de Beta (et donc phi)
		for (int r=0; r<k; r++)
		{
			//Compute Xr = X(Z==r,:) and Yr = Y(Z==r,:)
			int cardClustR=0;
			for (int i=0; i<n; i++)
			{
				if (Z[i] == r)
				{
					for (int j=0; j<p; j++)
						Xr[mi(cardClustR,j,n,p)] = X[mi(i,j,n,p)];
					for (int j=0; j<m; j++)
						Yr[mi(cardClustR,j,n,m)] = Y[mi(i,j,n,m)];
					cardClustR++;
				}
			}
			if (cardClustR == 0)
				continue;

			//Compute tXrXr = t(Xr) * Xr
			for (int j=0; j<p; j++)
			{
				for (int jj=0; jj<p; jj++)
				{
					double 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;
				}
			}

			//Get pseudo inverse = (t(Xr)*Xr)^{-1}
			double* 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;
					for (int u=0; u<cardClustR; u++)
						dotProduct += Xr[mi(u,j,n,p)] * Yr[mi(u,j,n,m)];
					tXrYr[mi(j,jj,p,m)] = dotProduct;
				}
			}

			//Fill matrixM with inverse * tXrYr = (t(Xr)*Xr)^{-1} * t(Xr) * Yr
			for (int j=0; j<p; j++)
			{
				for (int jj=0; jj<m; jj++)
				{
					double 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;
				}
			}
			free(invTXrXr);

			//U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr
			gsl_linalg_SV_decomp(matrixM, V, S, work);

			//Set m-rank(r) singular values to zero, and recompose
			//best rank(r) approximation of the initial product
			for (int j=rank[r]; j<m; j++)
				S->data[j] = 0.0;
			
			//[intermediate step] Compute hatBetaR = U * S * t(V)
			double* U = matrixM->data;
			for (int j=0; j<p; j++)
			{
				for (int jj=0; jj<m; jj++)
				{
					double 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;
				}
			}

			//Compute phi(:,:,r) = hatBetaR * Rho(:,:,r)
			for (int j=0; j<p; j++)
			{
				for (int jj=0; jj<m; jj++)
				{
					double 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;
				}
		}
		}

		/////////////
		// Etape E //
		/////////////
		
		double sumLogLLF2 = 0.0;
		for (int i=0; i<n; i++)
		{
			double sumLLF1 = 0.0;
			double maxLogGamIR = -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,:)...
				//*phi(:,:,r)) * transpose( Y(i,:)*Rho(:,:,r) - X(i,:)*phi(:,:,r) ) );
				//split in several sub-steps
				
				//compute det(Rho(:,:,r)) [TODO: avoid re-computations]
				for (int j=0; j<m; j++)
				{
					for (int jj=0; jj<m; jj++)
						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);

				//compute Y(i,:)*Rho(:,:,r)
				for (int j=0; j<m; j++)
				{
					YiRhoR[j] = 0.0;
					for (int u=0; u<m; u++)
						YiRhoR[j] += Y[mi(i,u,n,m)] * Rho[ai(u,j,r,m,m,k)];
				}

				//compute X(i,:)*phi(:,:,r)
				for (int j=0; j<m; j++)
				{
					XiPhiR[j] = 0.0;
					for (int u=0; u<p; u++)
						XiPhiR[j] += X[mi(i,u,n,p)] * phi[ai(u,j,r,p,m,k)];
				}

				//compute dotProduct < Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) . Y(:,i)*rho(:,:,r)-X(i,:)*phi(:,:,r) >
				double 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;

				//Z(i) = index of max (gam(i,:))
				if (logGamIR > maxLogGamIR)
				{
					Z[i] = r;
					maxLogGamIR = logGamIR;
				}
				sumLLF1 += exp(logGamIR) / pow(2*M_PI,m/2.0);
			}

			sumLogLLF2 += log(sumLLF1);
		}

		// Assign output variable LLF
		*LLF = -invN * sumLogLLF2;

		//newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi))));
		double 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)])
						/ (1.0+fabs(phi[ai(j,jj,r,p,m,k)]));
					if (tmpDist > newDeltaPhi)
						newDeltaPhi = tmpDist;
				}
			}
		}

		//update distance parameter to check algorithm convergence (delta(phi, Phi))
		//TODO: deltaPhi should be a linked list for perf.
		if (ite < deltaPhiBufferSize)
			deltaPhi[ite] = newDeltaPhi;
		else
		{
			sumDeltaPhi -= deltaPhi[0];
			for (int u=0; u<deltaPhiBufferSize-1; u++)
				deltaPhi[u] = deltaPhi[u+1];
			deltaPhi[deltaPhiBufferSize-1] = newDeltaPhi;
		}
		sumDeltaPhi += newDeltaPhi;

		// update other local variables
		for (int j=0; j<m; j++)
		{
			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)];
			}
		}
		ite++;
	}

	//free memory
	free(hatBetaR);
	free(deltaPhi);
	free(Phi);
	gsl_matrix_free(matrixE);
	gsl_matrix_free(matrixM);
	gsl_permutation_free(permutation);
	gsl_vector_free(work);
	gsl_matrix_free(V);
	gsl_vector_free(S);
	free(XiPhiR);
	free(YiRhoR);
	free(Xr);
	free(Yr);
	free(tXrXr);
	free(tXrYr);
	free(Z);
}
Commit	Line	Data
1d3c1faa BA	1	#include "EMGrank.h"
	2	#include <gsl/gsl_linalg.h>
	3
	4	// Compute pseudo-inverse of a square matrix
552b00e2	5	static double* pinv(const double* matrix, int dim)
1d3c1faa BA	6	{
	7	gsl_matrix* U = gsl_matrix_alloc(dim,dim);
	8	gsl_matrix* V = gsl_matrix_alloc(dim,dim);
	9	gsl_vector* S = gsl_vector_alloc(dim);
	10	gsl_vector* work = gsl_vector_alloc(dim);
552b00e2	11	double EPS = 1e-10; //threshold for singular value "== 0"
1d3c1faa BA	12
1d3c1faa BA	13	//copy matrix into U
552b00e2 BA	14	copyArray(matrix, U->data, dim*dim);
552b00e2 BA	15
1d3c1faa BA	16	//U,S,V = SVD of matrix
	17	gsl_linalg_SV_decomp(U, V, S, work);
	18	gsl_vector_free(work);
552b00e2	19
1d3c1faa	20	// Obtain pseudo-inverse by VS^{-1}t(U)
552b00e2 BA	21	double* inverse = (double)malloc(dimdim*sizeof(double));
552b00e2 BA	22	for (int i=0; i<dim; i++)
1d3c1faa	23	{
552b00e2	24	for (int ii=0; ii<dim; ii++)
1d3c1faa	25	{
552b00e2 BA	26	double dotProduct = 0.0;
552b00e2 BA	27	for (int j=0; j<dim; j++)
1d3c1faa BA	28	dotProduct += V->data[idim+j] (S->data[j] > EPS ? 1.0/S->data[j] : 0.0) * U->data[ii*dim+j];
	29	inverse[i*dim+ii] = dotProduct;
	30	}
	31	}
	32
	33	gsl_matrix_free(U);
	34	gsl_matrix_free(V);
	35	gsl_vector_free(S);
	36	return inverse;
	37	}
	38
	39	// TODO: comment EMGrank purpose
	40	void EMGrank(
	41	// IN parameters
552b00e2 BA	42	const double* Pi, // parametre de proportion
	43	const double* Rho, // parametre initial de variance renormalisé
	44	int mini, // nombre minimal d'itérations dans l'algorithme EM
	45	int maxi, // nombre maximal d'itérations dans l'algorithme EM
	46	const double* X, // régresseurs
	47	const double* Y, // réponse
	48	double tau, // seuil pour accepter la convergence
	49	const int* rank, // vecteur des rangs possibles
1d3c1faa	50	// OUT parameters
552b00e2 BA	51	double* phi, // parametre de moyenne renormalisé, calculé par l'EM
552b00e2 BA	52	double* LLF, // log vraisemblance associé à cet échantillon, pour les valeurs estimées des paramètres
1d3c1faa	53	// additional size parameters
552b00e2 BA	54	int n, // taille de l'echantillon
	55	int p, // nombre de covariables
	56	int m, // taille de Y (multivarié)
	57	int k) // nombre de composantes
1d3c1faa BA	58	{
1d3c1faa BA	59	// Allocations, initializations
552b00e2 BA	60	double* Phi = (double)calloc(pm*k,sizeof(double));
552b00e2 BA	61	double* hatBetaR = (double)malloc(pm*sizeof(double));
1d3c1faa	62	int signum;
552b00e2	63	double invN = 1.0/n;
1d3c1faa	64	int deltaPhiBufferSize = 20;
552b00e2 BA	65	double* deltaPhi = (double)malloc(deltaPhiBufferSizesizeof(double));
	66	int ite = 0;
	67	double sumDeltaPhi = 0.0;
	68	double* YiRhoR = (double)malloc(msizeof(double));
	69	double* XiPhiR = (double)malloc(msizeof(double));
	70	double* Xr = (double)malloc(np*sizeof(double));
	71	double* Yr = (double)malloc(nm*sizeof(double));
	72	double* tXrXr = (double)malloc(pp*sizeof(double));
	73	double* tXrYr = (double)malloc(pm*sizeof(double));
	74	gsl_matrix* matrixM = gsl_matrix_alloc(p, m);
1d3c1faa BA	75	gsl_matrix* matrixE = gsl_matrix_alloc(m, m);
	76	gsl_permutation* permutation = gsl_permutation_alloc(m);
	77	gsl_matrix* V = gsl_matrix_alloc(m,m);
	78	gsl_vector* S = gsl_vector_alloc(m);
	79	gsl_vector* work = gsl_vector_alloc(m);
	80
	81	//Initialize class memberships (all elements in class 0; TODO: randomize ?)
552b00e2	82	int* Z = (int*)calloc(n, sizeof(int));
1d3c1faa BA	83
1d3c1faa BA	84	//Initialize phi to zero, because some M loops might exit before phi affectation
552b00e2	85	for (int i=0; i<pmk; i++)
1d3c1faa BA	86	phi[i] = 0.0;
	87
	88	while (ite<mini \|\| (ite<maxi && sumDeltaPhi>tau))
552b00e2	89	{
1d3c1faa BA	90	/////////////
	91	// Etape M //
	92	/////////////
	93
	94	//M step: Mise à jour de Beta (et donc phi)
552b00e2	95	for (int r=0; r<k; r++)
1d3c1faa BA	96	{
1d3c1faa BA	97	//Compute Xr = X(Z==r,:) and Yr = Y(Z==r,:)
552b00e2 BA	98	int cardClustR=0;
552b00e2 BA	99	for (int i=0; i<n; i++)
1d3c1faa BA	100	{
	101	if (Z[i] == r)
	102	{
552b00e2 BA	103	for (int j=0; j<p; j++)
	104	Xr[mi(cardClustR,j,n,p)] = X[mi(i,j,n,p)];
	105	for (int j=0; j<m; j++)
	106	Yr[mi(cardClustR,j,n,m)] = Y[mi(i,j,n,m)];
1d3c1faa BA	107	cardClustR++;
	108	}
	109	}
552b00e2	110	if (cardClustR == 0)
1d3c1faa BA	111	continue;
	112
	113	//Compute tXrXr = t(Xr) * Xr
552b00e2	114	for (int j=0; j<p; j++)
1d3c1faa	115	{
552b00e2	116	for (int jj=0; jj<p; jj++)
1d3c1faa	117	{
552b00e2 BA	118	double dotProduct = 0.0;
	119	for (int u=0; u<cardClustR; u++)
	120	dotProduct += Xr[mi(u,j,n,p)] * Xr[mi(u,jj,n,p)];
	121	tXrXr[mi(j,jj,p,p)] = dotProduct;
1d3c1faa BA	122	}
	123	}
	124
	125	//Get pseudo inverse = (t(Xr)*Xr)^{-1}
552b00e2	126	double* invTXrXr = pinv(tXrXr, p);
1d3c1faa BA	127
1d3c1faa BA	128	// Compute tXrYr = t(Xr) * Yr
552b00e2	129	for (int j=0; j<p; j++)
1d3c1faa	130	{
552b00e2	131	for (int jj=0; jj<m; jj++)
1d3c1faa	132	{
552b00e2 BA	133	double dotProduct = 0.0;
	134	for (int u=0; u<cardClustR; u++)
	135	dotProduct += Xr[mi(u,j,n,p)] * Yr[mi(u,j,n,m)];
	136	tXrYr[mi(j,jj,p,m)] = dotProduct;
1d3c1faa BA	137	}
	138	}
	139
	140	//Fill matrixM with inverse * tXrYr = (t(Xr)Xr)^{-1} t(Xr) * Yr
552b00e2	141	for (int j=0; j<p; j++)
1d3c1faa	142	{
552b00e2	143	for (int jj=0; jj<m; jj++)
1d3c1faa	144	{
552b00e2 BA	145	double dotProduct = 0.0;
	146	for (int u=0; u<p; u++)
	147	dotProduct += invTXrXr[mi(j,u,p,p)] * tXrYr[mi(u,jj,p,m)];
	148	matrixM->data[j*m+jj] = dotProduct;
1d3c1faa BA	149	}
	150	}
	151	free(invTXrXr);
552b00e2	152
1d3c1faa BA	153	//U,S,V = SVD of (t(Xr)Xr)^{-1} * t(Xr) * Yr
1d3c1faa BA	154	gsl_linalg_SV_decomp(matrixM, V, S, work);
552b00e2 BA	155
552b00e2 BA	156	//Set m-rank(r) singular values to zero, and recompose
1d3c1faa	157	//best rank(r) approximation of the initial product
552b00e2	158	for (int j=rank[r]; j<m; j++)
1d3c1faa BA	159	S->data[j] = 0.0;
	160
	161	//[intermediate step] Compute hatBetaR = U * S * t(V)
552b00e2 BA	162	double* U = matrixM->data;
552b00e2 BA	163	for (int j=0; j<p; j++)
1d3c1faa	164	{
552b00e2	165	for (int jj=0; jj<m; jj++)
1d3c1faa	166	{
552b00e2 BA	167	double dotProduct = 0.0;
552b00e2 BA	168	for (int u=0; u<m; u++)
1d3c1faa	169	dotProduct += U[jm+u] S->data[u] * V->data[jj*m+u];
552b00e2	170	hatBetaR[mi(j,jj,p,m)] = dotProduct;
1d3c1faa BA	171	}
1d3c1faa BA	172	}
552b00e2	173
1d3c1faa	174	//Compute phi(:,:,r) = hatBetaR * Rho(:,:,r)
552b00e2	175	for (int j=0; j<p; j++)
1d3c1faa	176	{
552b00e2	177	for (int jj=0; jj<m; jj++)
1d3c1faa	178	{
552b00e2 BA	179	double dotProduct=0.0;
	180	for (int u=0; u<m; u++)
	181	dotProduct += hatBetaR[mi(j,u,p,m)] * Rho[ai(u,jj,r,m,m,k)];
	182	phi[ai(j,jj,r,p,m,k)] = dotProduct;
1d3c1faa	183	}
552b00e2	184	}
1d3c1faa BA	185	}
	186
	187	/////////////
	188	// Etape E //
	189	/////////////
	190
552b00e2 BA	191	double sumLogLLF2 = 0.0;
552b00e2 BA	192	for (int i=0; i<n; i++)
1d3c1faa	193	{
552b00e2 BA	194	double sumLLF1 = 0.0;
	195	double maxLogGamIR = -INFINITY;
	196	for (int r=0; r<k; r++)
1d3c1faa BA	197	{
	198	//Compute
	199	//Gam(i,r) = Pi(r) * det(Rho(:,:,r)) * exp( -1/2 * (Y(i,:)*Rho(:,:,r) - X(i,:)...
552b00e2	200	//phi(:,:,r)) transpose( Y(i,:)Rho(:,:,r) - X(i,:)phi(:,:,r) ) );
1d3c1faa BA	201	//split in several sub-steps
	202
	203	//compute det(Rho(:,:,r)) [TODO: avoid re-computations]
552b00e2	204	for (int j=0; j<m; j++)
1d3c1faa	205	{
552b00e2 BA	206	for (int jj=0; jj<m; jj++)
552b00e2 BA	207	matrixE->data[j*m+jj] = Rho[ai(j,jj,r,m,m,k)];
1d3c1faa BA	208	}
1d3c1faa BA	209	gsl_linalg_LU_decomp(matrixE, permutation, &signum);
552b00e2 BA	210	double detRhoR = gsl_linalg_LU_det(matrixE, signum);
552b00e2 BA	211
1d3c1faa	212	//compute Y(i,:)*Rho(:,:,r)
552b00e2	213	for (int j=0; j<m; j++)
1d3c1faa BA	214	{
1d3c1faa BA	215	YiRhoR[j] = 0.0;
552b00e2 BA	216	for (int u=0; u<m; u++)
552b00e2 BA	217	YiRhoR[j] += Y[mi(i,u,n,m)] * Rho[ai(u,j,r,m,m,k)];
1d3c1faa	218	}
552b00e2	219
1d3c1faa	220	//compute X(i,:)*phi(:,:,r)
552b00e2	221	for (int j=0; j<m; j++)
1d3c1faa BA	222	{
1d3c1faa BA	223	XiPhiR[j] = 0.0;
552b00e2 BA	224	for (int u=0; u<p; u++)
552b00e2 BA	225	XiPhiR[j] += X[mi(i,u,n,p)] * phi[ai(u,j,r,p,m,k)];
1d3c1faa	226	}
552b00e2	227
1d3c1faa	228	//compute dotProduct < Y(:,i)rho(:,:,r)-X(i,:)phi(:,:,r) . Y(:,i)rho(:,:,r)-X(i,:)phi(:,:,r) >
552b00e2 BA	229	double dotProduct = 0.0;
552b00e2 BA	230	for (int u=0; u<m; u++)
1d3c1faa	231	dotProduct += (YiRhoR[u]-XiPhiR[u]) * (YiRhoR[u]-XiPhiR[u]);
552b00e2 BA	232	double logGamIR = log(Pi[r]) + log(detRhoR) - 0.5*dotProduct;
552b00e2 BA	233
1d3c1faa BA	234	//Z(i) = index of max (gam(i,:))
	235	if (logGamIR > maxLogGamIR)
	236	{
	237	Z[i] = r;
	238	maxLogGamIR = logGamIR;
	239	}
	240	sumLLF1 += exp(logGamIR) / pow(2*M_PI,m/2.0);
	241	}
552b00e2	242
1d3c1faa BA	243	sumLogLLF2 += log(sumLLF1);
1d3c1faa BA	244	}
552b00e2	245
1d3c1faa BA	246	// Assign output variable LLF
	247	LLF = -invN sumLogLLF2;
	248
	249	//newDeltaPhi = max(max((abs(phi-Phi))./(1+abs(phi))));
552b00e2 BA	250	double newDeltaPhi = 0.0;
552b00e2 BA	251	for (int j=0; j<p; j++)
1d3c1faa	252	{
552b00e2	253	for (int jj=0; jj<m; jj++)
1d3c1faa	254	{
552b00e2	255	for (int r=0; r<k; r++)
1d3c1faa	256	{
552b00e2 BA	257	double tmpDist = fabs(phi[ai(j,jj,r,p,m,k)]-Phi[ai(j,jj,r,p,m,k)])
552b00e2 BA	258	/ (1.0+fabs(phi[ai(j,jj,r,p,m,k)]));
1d3c1faa BA	259	if (tmpDist > newDeltaPhi)
	260	newDeltaPhi = tmpDist;
	261	}
	262	}
	263	}
	264
	265	//update distance parameter to check algorithm convergence (delta(phi, Phi))
	266	//TODO: deltaPhi should be a linked list for perf.
	267	if (ite < deltaPhiBufferSize)
	268	deltaPhi[ite] = newDeltaPhi;
	269	else
	270	{
	271	sumDeltaPhi -= deltaPhi[0];
	272	for (int u=0; u<deltaPhiBufferSize-1; u++)
	273	deltaPhi[u] = deltaPhi[u+1];
	274	deltaPhi[deltaPhiBufferSize-1] = newDeltaPhi;
	275	}
	276	sumDeltaPhi += newDeltaPhi;
	277
	278	// update other local variables
552b00e2	279	for (int j=0; j<m; j++)
1d3c1faa	280	{
552b00e2	281	for (int jj=0; jj<p; jj++)
1d3c1faa	282	{
552b00e2 BA	283	for (int r=0; r<k; r++)
552b00e2 BA	284	Phi[ai(j,jj,r,p,m,k)] = phi[ai(j,jj,r,p,m,k)];
1d3c1faa BA	285	}
1d3c1faa BA	286	}
552b00e2	287	ite++;
1d3c1faa	288	}
552b00e2	289
1d3c1faa BA	290	//free memory
	291	free(hatBetaR);
	292	free(deltaPhi);
	293	free(Phi);
	294	gsl_matrix_free(matrixE);
	295	gsl_matrix_free(matrixM);
	296	gsl_permutation_free(permutation);
	297	gsl_vector_free(work);
	298	gsl_matrix_free(V);
	299	gsl_vector_free(S);
	300	free(XiPhiR);
	301	free(YiRhoR);
	302	free(Xr);
	303	free(Yr);
	304	free(tXrXr);
	305	free(tXrYr);
552b00e2	306	free(Z);
1d3c1faa	307	}