// PLSAT.cpp : Defines the entry point for the console application.
// Quick DEMO of non-negative matrix factorization and 
// probabilistic latent semantic analysis. Andrew Polar, ezcodesample.com, semanticsearchart.com.
// Oct. 2011

#include "stdafx.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <memory.h>

const int numberOfCategories = 2; 
const int numberOfRows       = 5;  
const int numberOfColumns    = 6;  
static float data[numberOfRows][numberOfColumns] = {
	{9.0, 2.0, 1.0, 1.0, 1.0, 0.0},   
	{8.0, 3.0, 2.0, 1.0, 0.0, 0.0},   
	{3.0, 0.0, 0.0, 1.0, 2.0, 8.0},
	{0.0, 1.0, 0.0, 2.0, 4.0, 7.0},
	{2.0, 1.0, 1.0, 0.0, 1.0, 3.0},   
};

//////////////Reusable functions/////////////////////
void printVector(char* title, float* Vector, int n) {
	printf("%s\n\n", title);
	for (int i=0; i<n; ++i) {
		printf("%10.4f", Vector[i]);
	}
	printf("\n");
}

void printMatrix(char* title, float** Matrix, int nRows, int nCols) {
	printf("%s\n\n", title);
	for (int i=0; i<nRows; ++i) {
		for (int j=0; j<nCols; ++j) {
			printf("%10.4f", Matrix[i][j]);
		}
		printf("\n");
	}
	printf("\n");
}

void normalizeMatrix(float** Matrix, int nRows, int nCols) {
	for (int i=0; i<nRows; ++i) {
		float min = Matrix[i][0];
		for (int j=0; j<nCols; ++j) {
			if (min > Matrix[i][j]) min = Matrix[i][j];
		}
		if (min < 0.0) {
			for (int j=0; j<nCols; ++j) {
				Matrix[i][j] -= min;
			}
		}
		float s = 0.0;
		for (int j=0; j<nCols; ++j) {
			s += Matrix[i][j];
		}
		if (s > 0.0) {
			for (int j=0; j<nCols; ++j) {
				Matrix[i][j] /= s;
			}
		}
	}
}

void initializeMatrix(float** Matrix, int nRows, int nCols) {
	srand((unsigned)(time(0))); 
	for (int i=0; i<nRows; ++i) {
		for (int j=0; j<nCols; ++j) {
			Matrix[i][j] = (float)(rand() % 100);
			Matrix[i][j] /= 100.0;
		}
	}
}

void initializeZ(float* Z, int nSize) {
	for (int i=0; i<nSize; ++i) Z[i] = 1.0/(float)(nSize);
}

float computeLikelihood(float** D, float** W, int nRows, int nCols, int nCats) {
	float likelihood = 0.0;
	for (int i=0; i<nRows; ++i) {
		for (int j=0; j<nCols; ++j) {
			if (data[i][j] > 0.0) {
				float s = 0.0;
				for (int k = 0; k<nCats; ++k) {
					s += D[i][k] * W[k][j];
				}
				if (s > 0.0) {
					likelihood += data[i][j] * log(s);
				}
			}
		}
	}
	return likelihood;
}
////////End///////////////////////////


//Next is non-negative matrix factorization
//Cholesky matrix inversion
void choldc1(int n, float** a, float* p) {
	int i,j,k;
    float sum;
	for (i = 0; i<n; i++) {
		for (j = i; j<n; j++) {
			sum = a[i][j];
            for (k = i - 1; k >= 0; k--) {
				sum -= a[i][k] * a[j][k];
			}
            if (i == j) {
				if(sum <= 0) { 
					p[i] = 1.0; //regularization
				}
				else {
					p[i] = sqrt(sum);
				}
			}
            else {
				a[j][i] = sum / p[i];
			}
		}
	}
}

void choldc(int n, float** A, float** a) {
	int i,j;
	float* p = (float*)malloc(n * sizeof(float));
    for (i = 0; i<n; i++) 
		for (j = 0; j<n; j++) 
			a[i][j] = A[i][j];
	choldc1(n, a, p);
	for (i = 0; i<n; i++) {
		a[i][i] = p[i];
        for (j = i + 1; j<n; j++) {
			a[i][j] = 0;
		}
	}
	if(p) free(p);
}

//inverts positive definite M, returns result in M 
void choleskyInvert(float** M, int size) {
	float** T = (float**)malloc(size * sizeof(float*));
	for(int i = 0; i<size; ++i) {
		T[i] = (float*)malloc(size * sizeof(float));
	}

	choldc(size, M, T); //find lower triangular
	
	//make it symmetric
	for (int i = 0; i<size; ++i)
		for (int j = i+1; j<size; ++j)
			T[i][j] = T[j][i];

	float* v = (float*)malloc(size * sizeof(float));
	float* e = (float*)malloc(size * sizeof(float));
	
	for (int loop=0; loop<size; ++loop) {
		memset(e, 0x00, size * sizeof(float));
		e[loop] = 1.0;

		//solve lower triangular
		for (int i = 0; i<size; ++i) {
			float sum = 0.0;
			for (int j = 0; j<i; ++j) {
				sum += T[i][j] * v[j];
			}
			v[i] = (e[i] - sum) / T[i][i];
		}

		//solve upper triangular
		for (int i = size - 1; i>= 0; --i) {
			float sum = 0.0;
			for (int j = size - 1; j>i; --j) {
				sum += T[i][j] * M[j][loop];
			}
			M[i][loop] = (v[i] - sum) / T[i][i];
		}
	}

	if (e) free(e);
	if (v) free(v);
	if (T) {
		for (int i = 0; i<size; ++i) {
			if (T[i])free(T[i]);
		}
		free(T);
	}
}
//end Cholesky inversion

void makeN(float** N, int nRows, int nCols) {
	for (int i=0; i<nRows; ++i) {
		for (int j=0; j<nCols; ++j) {
			N[i][j] = data[i][j];
		}
	}
	normalizeMatrix(N, nRows, nCols);
}

void getNewW(float** D, float** W, float** N, int nRows, int nCols, int nCats) {
	//it is D = N * WT multiplication
	for (int nWhich=0; nWhich<nCats; ++nWhich) {
		for (int i=0; i<nRows; ++i) {
			D[i][nWhich] = 0.0;
			for (int j=0; j<nCols; ++j) {
				if (data[i][j] > 0.0) {
					D[i][nWhich] += N[i][j] * W[nWhich][j]; 
				}
			}
		}
	}
	//copy D into N
	for (int i=0; i<nRows; ++i) {
		for (int j=0; j<nCats; ++j) {
			N[i][j] = D[i][j];
		}
	}
	//it is WWT computation
	float** WWT = (float**)malloc(nCats * sizeof(float*));
	for (int i=0; i<nCats; ++i) {
		WWT[i] = (float*)malloc(nCats * sizeof(float));
	}
	for (int i=0; i<nCats; ++i) {
		for (int j=0; j<nCats; ++j) {
			WWT[i][j] = 0.0;
			for (int k=0; k<nCols; ++k) {
				WWT[i][j] += W[i][k] * W[j][k];
			}
		}
	}
	choleskyInvert(WWT, nCats);
	//it is inverse(WWT) * D multiplication
	for (int i=0; i<nRows; ++i) {
		for (int j=0; j<nCats; ++j) {
			D[i][j] = 0.0;
			for (int k=0; k<nCats; ++k) {
				D[i][j] += N[i][k] * WWT[k][j];
			}
		}
	}
	normalizeMatrix(D, nRows, nCats);
	if (WWT) {
		for (int i=0; i<nCats; ++i) {
			if (WWT[i]) free(WWT[i]);
		}
		free(WWT);
	}
}

void getNewD(float** D, float** W, float** N, int nRows, int nCols, int nCats) {
	//it is W = DT * N multiplication
	for(int nWhich=0; nWhich<nCats; ++nWhich) {
		for (int j=0; j<nCols; ++j) {
			W[nWhich][j] = 0.0;
			for (int i=0; i<nRows; ++i) {
				if (N[i][j] > 0.0) {
					W[nWhich][j] += N[i][j] * D[i][nWhich]; 
				}
			}
		}
	}
	//copy W into N
	for (int i=0; i<nCats; ++i) {
		for (int j=0; j<nCols; ++j) {
			N[i][j] = W[i][j];
		}
	}
	//it is DTD computation
	float** DTD = (float**)malloc(nCats * sizeof(float*));
	for (int i=0; i<nCats; ++i) {
		DTD[i] = (float*)malloc(nCats * sizeof(float));
	}
	for (int i=0; i<nCats; ++i) {
		for (int j=0; j<nCats; ++j) {
			DTD[i][j] = 0.0;
			for (int k=0; k<nRows; ++k) {
				DTD[i][j] += D[k][i] * D[k][j];
			}
		}
	}
	choleskyInvert(DTD, nCats);
	//it is inverse(DTD) * W multiplication
	for (int i=0; i<nCols; ++i) {
		for (int j=0; j<nCats; ++j) {
			W[j][i] = 0.0;
			for (int k=0; k<nCats; ++k) {
				W[j][i] += DTD[k][j] * N[k][i];
			}
		}
	}
	normalizeMatrix(W, nCats, nCols);
	if (DTD) {
		for (int i=0; i<nCats; ++i) {
			if (DTD[i]) free(DTD[i]);
		}
		free(DTD);
	}
}

bool makeApproximationStepNMF(float** D, float** W, float** N, int nRows, int nCols, int nCats) {
	static float prevLikelihood = 0.0;
	makeN(N, nRows, nCols);
	getNewD(D, W, N, nRows, nCols, nCats);
	makeN(N, nRows, nCols);
	getNewW(D, W, N, nRows, nCols, nCats);
	float likelihood = computeLikelihood(D, W, nRows, nCols, nCats);
	printf("current likelihood %f\n", likelihood);
	if (fabs(prevLikelihood  - likelihood) / fabs(likelihood) < 0.001) {
		return true;
	}
	prevLikelihood = likelihood;
	return false;
}

void showNonNegativeMatrixFactoring() {
	int nRows = numberOfRows;
	int nCols = numberOfColumns;
	int nCats = numberOfCategories;

	//Allocate memory
	float** D = (float**)malloc(nRows * sizeof(float*));
	for (int i=0; i<nRows; ++i) {
		D[i] = (float*)malloc(nCats * sizeof(float));
	}

	float** W = (float**)malloc(nCats * sizeof(float*));
	for (int i=0; i<nCats; ++i) {
		W[i] = (float*)malloc(nCols * sizeof(float));
	}

	float** N = (float**)malloc(nRows * sizeof(float*));
	for (int i=0; i<nRows; ++i) {
		N[i] = (float*)malloc(nCols * sizeof(float));
	}
	//end

	/////////////////////////////Start pLSA/////////////////////////////
	//compute ideal likelihood
	makeN(N, nRows, nCols);
	float theoreticalLimitOfLikelihood = 0.0;
	for (int i=0; i<nRows; ++i) {
		for (int j=0; j<nCols; ++j) {
			if (data[i][j] > 0) {
				theoreticalLimitOfLikelihood += data[i][j] * log(N[i][j]);
			}
		}
	}
	printf("Start pLSA NMF, theoretical limit for likelihood %f\n\n", theoreticalLimitOfLikelihood);

	initializeMatrix(D, nRows, nCats);
	initializeMatrix(W, nCats, nCols);
	normalizeMatrix(D, nRows, nCats);
	normalizeMatrix(W, nCats, nCols);
	while (!makeApproximationStepNMF(D, W, N, nRows, nCols, nCats)) {}
	printf("\n");
	printMatrix("Document-category", D, nRows, nCats);
	printMatrix("Category-word", W, nCats, nCols);

	printf("End pLSA NMF \n\n");
	/////////////////////////////End pLSA///////////////////////////////

	//Deallocate memory
	if (N) {
		for (int i=0; i<nRows; ++i) {
			if (N[i]) free(N[i]);
		}
		free(N);
	}

	if (W) {
		for (int i=0; i<nCats; i++) {
			if (W[i]) free(W[i]);
		}
		free(W);
	}
	if (D) {
		for (int i=0; i<nRows; i++) {
			if (D[i]) free(D[i]);
		}
		free(D);
	}
	//end
}
//End non-negative matrix factorization

//Next is PLSA
void getNewN(float** N, float** D, float** W, float* Z, int nRows, int nCols, int nCats) {
	for (int i=0; i<nRows; ++i) {
		for (int j=0; j<nCols; ++j) {
			if (data[i][j] > 0.0) {
				float s = 0.0;
				for (int k = 0; k<nCats; ++k) {
					s += D[i][k] * W[k][j] * Z[k];
				}
				if (s > 0.0) {
					N[i][j] = data[i][j] / s;
				}
				else {
					N[i][j] = 0.0;
				}
			}
			else {
				N[i][j] = 0.0;
			}
		}
	}
}

void getNewD(float** D2, float** D1, float** W2, float** W1, float** N, float* Z, int nRows, int nCols, int nCats) {
	//it is D = N * WT multiplication
	for (int nWhich=0; nWhich<nCats; ++nWhich) {
		for (int i=0; i<nRows; ++i) {
			D2[i][nWhich] = 0.0;
			for (int j=0; j<nCols; ++j) {
				if (data[i][j] > 0.0) {
					D2[i][nWhich] += N[i][j] * W1[nWhich][j]; 
				}
			}
		}
	}
	for (int i=0; i<nRows; ++i) {
		for (int nWhich=0; nWhich<nCats; ++nWhich) {
			D2[i][nWhich] *= D1[i][nWhich];
		}
	}
}

void getNewW(float** W2, float** W1, float** D2, float** D1, float** N, float* Z, int nRows, int nCols, int nCats) {
	//it is W = DT * N multiplication
	for(int nWhich=0; nWhich<nCats; ++nWhich) {
		for (int j=0; j<nCols; ++j) {
			W2[nWhich][j] = 0.0;
			for (int i=0; i<nRows; ++i) {
				if (data[i][j] > 0.0) {
					W2[nWhich][j] += N[i][j] * D1[i][nWhich]; 
				}
			}
		}
	}
	for(int nWhich=0; nWhich<nCats; ++nWhich) {
		for (int j=0; j<nCols; ++j) {
			W2[nWhich][j] *= W1[nWhich][j]; 
		}
	}
}

void getNewZ(float* Z, float** D, int nRows, int nCols, int nCats) {
	for(int nWhich=0; nWhich<nCats; ++nWhich) {
		Z[nWhich] = 0.0;
		for (int j=0; j<nCols; ++j) {
			for (int i=0; i<nRows; ++i) {
				if (data[i][j] > 0.0) {
					Z[nWhich] += data[i][j] * D[i][nWhich]; 
				}
			}
		}
	}
	float s = 0.0;
	for (int i=0; i<nCats; ++i) {
		s += Z[i];
	}
	for (int i=0; i<nCats; ++i) {
		Z[i] /= s;
	}
}

bool makeApproximationStepPLSA(float** D1, float** D2, float** W1, float** W2, float** N, float* Z, int nRows, int nCols, int nCats) {
	//This is E-step
	getNewN(N, D1, W1, Z, nRows, nCols, nCats); 
	//Next is M-step
	getNewD(D2, D1, W2, W1, N, Z, nRows, nCols, nCats);
	getNewW(W2, W1, D2, D1, N, Z, nRows, nCols, nCats);
	getNewZ(Z, D1, nRows, nCols, nCats);
	
	normalizeMatrix(D2, nRows, nCats);
	normalizeMatrix(W2, nCats, nCols);

	float difference = 0.0;
	for (int i=0; i<nRows; ++i) {
		for (int j=0; j<nCats; ++j) {
			difference += fabs(D1[i][j] - D2[i][j]);
			D1[i][j] = D2[i][j];
		}
	}
	for (int i=0; i<nCats; ++i) {
		for (int j=0; j<nCols; ++j) {
			difference += fabs(W1[i][j] - W2[i][j]);
			W1[i][j] = W2[i][j];
		}
	}
	difference /= (float)(nCats);
	difference /= (float)(nCats);
	if (difference < 0.001) return true;
	return false;
}

void showPLSA() {

	int nRows = numberOfRows;
	int nCols = numberOfColumns;
	int nCats = numberOfCategories;

	//Allocate memory
	float** D1 = (float**)malloc(nRows * sizeof(float*));
	for (int i=0; i<nRows; ++i) {
		D1[i] = (float*)malloc(nCats * sizeof(float));
	}

	float** D2 = (float**)malloc(nRows * sizeof(float*));
	for (int i=0; i<nRows; ++i) {
		D2[i] = (float*)malloc(nCats * sizeof(float));
	}

	float** W1 = (float**)malloc(nCats * sizeof(float*));
	for (int i=0; i<nCats; ++i) {
		W1[i] = (float*)malloc(nCols * sizeof(float));
	}

	float** W2 = (float**)malloc(nCats * sizeof(float*));
	for (int i=0; i<nCats; ++i) {
		W2[i] = (float*)malloc(nCols * sizeof(float));
	}

	float** N = (float**)malloc(nRows * sizeof(float*));
	for (int i=0; i<nRows; ++i) {
		N[i] = (float*)malloc(nCols * sizeof(float));
	}

	float* Z = (float*)malloc(nCats * sizeof(float));
	//end

	/////////////////////////////Start pLSA/////////////////////////////
	printf("Start PLSA\n");
	initializeMatrix(D1, nRows, nCats);
	initializeMatrix(W1, nCats, nCols);
	normalizeMatrix(D1, nRows, nCats);
	normalizeMatrix(W1, nCats, nCols);
	initializeZ(Z, nCats);
	int nCounter = 0;
	while (!makeApproximationStepPLSA(D1, D2, W1, W2, N, Z, nRows, nCols, nCats)) {
		float currentLikelihood = computeLikelihood(D2, W2, nRows, nCols, nCats);
		printf("current likelihood %f\n", currentLikelihood);
		++nCounter;
		if (nCounter >= 300) break;
	}
	printf("\n\n");
	normalizeMatrix(D2, nRows, nCats);
	normalizeMatrix(W2, nCats, nCols);
	printMatrix("Document-category", D2, nRows, nCats);
	printMatrix("Word-category", W2, nCats, nCols);
	printf("End PLSA\n\n");
	/////////////////////////////End pLSA///////////////////////////////

	//Deallocate memory
	if (Z) free(Z);

	if (N) {
		for (int i=0; i<nRows; ++i) {
			if (N[i]) free(N[i]);
		}
		free(N);
	}

	if (W2) {
		for (int i=0; i<nCats; i++) {
			if (W2[i]) free(W2[i]);
		}
		free(W2);
	}
	if (W1) {
		for (int i=0; i<nCats; i++) {
			if (W1[i]) free(W1[i]);
		}
		free(W1);
	}
	if (D2) {
		for (int i=0; i<nRows; i++) {
			if (D2[i]) free(D2[i]);
		}
		free(D2);
	}
	if (D1) {
		for (int i=0; i<nRows; i++) {
			if (D1[i]) free(D1[i]);
		}
		free(D1);
	}
	//end
}
//End PLSA

int main() {

	showNonNegativeMatrixFactoring();
	printf("*************************************************\n\n");
	showPLSA();

	return 0;
}