NICE_VisLearning/math/cluster/GMM.cpp

#ifdef NICE_USELIB_OPENMP
#include <omp.h>
#endif

#include <stdio.h>

#include "GMM.h"
#include <core/vector/Algorithms.h>

#include "vislearning/math/cluster/KMeans.h"

#define DEBUGGMM

using namespace OBJREC;
using namespace NICE;
using namespace std;

GMM::GMM()
{
	gaussians = 3;
	maxiter = 200;
	featsperclass = 0;
	tau = 10.0;
	srand ( time (NULL) );
	comp = false;
}

GMM::GMM(int _gaussians):gaussians(_gaussians)
{
	maxiter = 200;
	featsperclass = 0;
	tau = 0.0;
	srand ( time (NULL) );
	pyramid = false;
	comp = false;
}

GMM::GMM(const Config *conf, int _gaussians):gaussians(_gaussians)
{
	this->conf = conf;
	
	if(gaussians < 2)
		gaussians = conf->gI("GMM", "gaussians", 2 );
		
	maxiter = conf->gI("GMM", "maxiter", 200 );
		
	featsperclass = conf->gI("GMM", "featsperclass", 0 );
	
	tau = conf->gD("GMM", "tau", 100.0);
	
	pyramid = conf->gB("GMM", "pyramid", false);
	
	comp = false;
	
	srand ( time (NULL) );
}

void GMM::computeMixture(Examples examples)
{
	int fsize = (int)examples.size();
	assert(fsize >= gaussians);
	
	dim = examples[0].second.vec->size();
	
	int samples = fsize;
	if(featsperclass > 0)
	{
		samples = featsperclass*gaussians;
		samples = std::min(samples, fsize);
	}
		
	// Copy data in Matrix
	VVector dataset;

	cout << "reduced training data for GMM from " << fsize << " features to " << samples << " features.";
	for(int i = 0; i < samples; i++)
	{
		int k = rand() % fsize;
		NICE::Vector *vec = examples[k].second.vec;
		dataset.push_back(*vec);
	}
	
	computeMixture(dataset);
}

void GMM::computeMixture(const VVector &DataSet)
{
	// Learn the GMM model
	assert(DataSet.size() >= (uint)gaussians);
	//initEMkMeans(DataSet); // initialize the model
	srand ( time (NULL) );
	
	bool poweroftwo = false;
	int power = 1;
	while(power <= gaussians)
	{
		if(gaussians == power)
			poweroftwo = true;
		power *= 2;
	}
	
	if(poweroftwo && pyramid)
	{
		initEM(DataSet); // initialize the model
		int g = 1;
		
		while(g <= gaussians)
		{
			cout << "g = " << g << endl;
			doEM(DataSet, g);
			
			if(2*g <= gaussians)
			{
				for(int i = g; i < g*2; i++)
				{
					mu[i] = mu[i-g];
					sparse_sigma[i] = sparse_sigma[i-g];
					sparse_inv_sigma[i] = sparse_inv_sigma[i-g];
					log_det_sigma[i] = log_det_sigma[i-g];
					priors[i] = priors[i-g];
					
					double interval = 1.0;
					for(int k = 0; k < (int)mu[i].size(); k++)
					{
						interval = mu[i][k];
						interval = std::max(interval, 0.1);
						double r = (interval*((double)rand()/(double)RAND_MAX))-interval/2.0;
						mu[i][k] += r;
					}
				}
			}
			g *= 2;
		}
	}
	else
	{
		initEMkMeans(DataSet); // initialize the model
		doEM(DataSet, gaussians);
	}
	// performs EM
}

inline double diagDeterminant(const NICE::Vector &sparse_mat)
{
	double det = 1.0;
	for(int i = 0; i < (int)sparse_mat.size(); i++)
	{
		det *= sparse_mat[i];
	}
	
	return det;
}

inline double logdiagDeterminant(const NICE::Vector &sparse_mat)
{
	double det = 0.0;
	for(int i = 0; i < (int)sparse_mat.size(); i++)
	{
		det += log(sparse_mat[i]);
	}

	return det;
}

inline NICE::Vector diagInverse(const NICE::Vector &sparse_mat)
{
	NICE::Vector inv(sparse_mat.size());
	for(int i = 0; i < (int)sparse_mat.size(); i++)
	{
		inv[i] = 1.0/sparse_mat[i];
	}
	return inv;
}

void GMM::initEMkMeans(const VVector &DataSet)
{
	/*init GMM with k-Means*/
	KMeans k( gaussians );
	VVector means;
	vector<double> weights;
	vector<int> assignment;
	k.cluster ( DataSet, means, weights, assignment );
	
	int nData = DataSet.size();
	this->dim = DataSet[0].size();
	cdimval = dim*log(2*3.14159);
	vector<int> pop(gaussians, 0);
	priors.resize(gaussians);
	mu = VVector(gaussians, dim);
	log_det_sigma.clear();
	vector<int> allk;
	
	NICE::Vector globmean(dim);
	globmean.set(0.0);
	
	for(int n = 0; n < (int)DataSet.size(); n++) /* getting the max value for time */
	{
		globmean = globmean + DataSet[n];
	}

	globmean *= (1.0/nData);
	
	NICE::Vector sparse_globsigma;

	sparse_globsigma.resize(dim);
	sparse_globsigma.set(0.0);
	
	for(int i = 0; i < (int)DataSet.size(); i++) // Covariances updates
	{
		// nur diagonal Elemente berechnen
		for(int d = 0; d < dim; d++)
		{
			double diff = (DataSet[i][d] - globmean[d]);
			sparse_globsigma[d] += diff*diff;
		}
	}
	
	sparse_globsigma *= (1.0/DataSet.size());
	
	for(int i = 0; i < gaussians; i++)
	{
		NICE::Vector _inv_sigma = diagInverse(sparse_globsigma);

		sparse_sigma.push_back(sparse_globsigma);
		sparse_inv_sigma.push_back(_inv_sigma);
		log_det_sigma.push_back(logdiagDeterminant(sparse_globsigma));
		mu[i] = means[i];
		//priors[i]=1.0/(double)gaussians; // set equi-probables states 
		priors[i]=weights[i]; // set kMeans weights 
	}
}

void GMM::initEM(const VVector &DataSet)
{

	/* init the GaussianMixture by using randomized meanvectors */
	int nData = DataSet.size();
	this->dim = DataSet[0].size();
	cdimval = dim*log(2*3.14159);
	vector<int> pop(gaussians, 0);
	priors.resize(gaussians);
	mu = VVector(gaussians, dim);
	log_det_sigma.clear();
	vector<int> allk;
	
	NICE::Vector globmean(dim);
	globmean.set(0.0);
	
	for(int n = 0; n < (int)DataSet.size(); n++) /* getting the max value for time */
	{
		globmean = globmean + DataSet[n];
	}

	globmean *= (1.0/nData);
	
	NICE::Vector sparse_globsigma;

	sparse_globsigma.resize(dim);
	sparse_globsigma.set(0.0);
	
	for(int i = 0; i < (int)DataSet.size(); i++) // Covariances updates
	{
		// nur diagonal Elemente berechnen
		for(int d = 0; d < dim; d++)
		{
			double diff = (DataSet[i][d] - globmean[d]);
			sparse_globsigma[d] += diff*diff;
		}
	}
	
	sparse_globsigma *= (1.0/DataSet.size());
	
	for(int i = 0; i < gaussians; i++)
	{
		priors[i]=1.0/(double)gaussians; // set equi-probables states 
		NICE::Vector _inv_sigma = diagInverse(sparse_globsigma);

		sparse_sigma.push_back(sparse_globsigma);
		sparse_inv_sigma.push_back(_inv_sigma);
		log_det_sigma.push_back(logdiagDeterminant(sparse_globsigma));
		bool newv = false;
		int k;
		while(!newv)
		{
			newv = true;
			k = rand() % nData;
		
			for(int nk = 0; nk < (int)allk.size(); nk++)
				if(allk[nk] == k)
				{
					newv = false;
				}
			if(newv)
				allk.push_back(k);
		}
		mu[i] = DataSet[k];
	}
}

inline void sumRow(NICE::Matrix mat, NICE::Vector &res)
{
	int row = mat.rows();
	int column = mat.cols();
	res =NICE::Vector(column);
	res.set(1e-5f);
//#pragma omp parallel for
	for(int i = 0; i < column; i++){
		for(int j = 0; j < row; j++){
			res[i] += mat(j, i);
		}
	}
}

double GMM::logpdfState(const NICE::Vector &Vin,int state)
{
	/* get the probability density for a given state and a given vector */
	double p;
	NICE::Vector dif;
		
	dif = Vin;
	dif -= mu[state];
		
	p = 0.0;
	for(int i = 0; i < (int)dif.size(); i++)
	{
		p += dif[i] * dif[i] * sparse_inv_sigma[state][i];
	}
	
	p = -0.5*(p + cdimval + log_det_sigma[state]);
	return p;
}

int GMM::doEM(const VVector &DataSet, int nbgaussians)
{
  /* perform Expectation/Maximization on the given Dataset :
	Matrix DataSet(nSamples,Dimensions).
	The GaussianMixture Object must be initialised before 
	(see initEM or initEMkMeans methods ) */
	
#ifdef DEBUG
	cerr << "GMM::start EM" << endl;
#endif	
	
	int nData = DataSet.size();
	int iter = 0;
	double log_lik;
	double log_lik_threshold = 1e-6f;
	double log_lik_old = -1e10f;

	NICE::Matrix unity(dim,dim,0.0);
	for(int k = 0; k < dim; k++) 
		unity(k, k)=1.0;
	
	//EM loop
	while(true)
	{
#ifdef DEBUGGMM
		cerr << "GMM::EM: iteration: " << iter << endl;
#endif
		
		vector<double> sum_p;
		sum_p.resize(nData);

		for(int i = 0; i < nData; i++)
		{
			sum_p[i] = 0.0;
		}
	
		NICE::Matrix pxi(nData,gaussians);
		pxi.set(0.0);
		NICE::Matrix pix(nData,gaussians);
		pix.set(0.0);
		NICE::Vector E;
		
		iter++;
		if (iter>maxiter){
			cerr << "EM stops here. Max number of iterations (" << maxiter << ") has been reached." << endl;
			return iter;
		}

		double sum_log = 0.0; 
		
    	// computing expectation
		double maxp = -numeric_limits<double>::max();
		vector<double> maxptmp(nData,-numeric_limits<double>::max());
#pragma omp parallel for
		for(int i = 0; i < nData; i++)
		{
			for(int j = 0; j < nbgaussians; j++)
			{
				double p = logpdfState(DataSet[i],j); // log(P(x|i))
				maxptmp[i] = std::max(maxptmp[i], p);
				pxi(i, j) = p;
			}
		}
		
		for(int i = 0; i < nData; i++)
		{
			maxp = std::max(maxp, maxptmp[i]);
		}

#pragma omp parallel for
		for(int i = 0; i < nData; i++)
		{
			sum_p[i] = 0.0;
			for(int j = 0; j < nbgaussians; j++)
			{
				double p = pxi(i, j) - maxp; // log(P(x|i))
				p = exp(p); // P(x|i)
				
				if(p < 1e-40)
					p = 1e-40;

				pxi(i, j) = p;
				sum_p[i] += p*priors[j];
			} 
		}
			
		for(int i = 0; i < nData; i++)
		{
			sum_log += log(sum_p[i]);
		}
		
#pragma omp parallel for
		for(int j = 0; j < nbgaussians; j++)
		{
			for(int i = 0; i < nData; i++)
			{
				pix(i, j) = (pxi(i, j)*priors[j])/sum_p[i]; // then P(i|x)			
			}   
		}

		// here we compute the log likehood 
		log_lik = sum_log/nData; 
#ifdef DEBUGGMM	
		cout << "diff: " << fabs((log_lik/log_lik_old)-1) << " thresh: " << log_lik_threshold << " sum: " << sum_log <<  endl;
		//cout << "logold: " << log_lik_old << " lognew: " << log_lik << endl;
#endif
		
		if(fabs((log_lik/log_lik_old)-1) < log_lik_threshold )
		{
			//if log likehood hasn't move enough, the algorithm has converged, exiting the loop 
			return iter;
		}

		log_lik_old = log_lik;

		// Update Step
		sumRow(pix,E);
		
#pragma omp parallel for
		for(int j = 0; j < nbgaussians; j++) 
		{
			priors[j] = (E[j]+tau)/(nData+tau*nbgaussians); // new priors 

			NICE::Vector tmu(dim);
			tmu.set(0.0);

			NICE::Vector sparse_tmsigma(dim);
			sparse_tmsigma.set(0.0);

			for(int i = 0; i < nData; i++) // Means update loop
			{
				tmu = tmu + (DataSet[i]*pix(i,j));
			}
									
			NICE::Vector tmu2 = mu[j];
			mu[j] = tmu + tau*tmu2;
			mu[j] = mu[j]*(1.0/(E[j]+tau));
			
			for(int i = 0; i < nData; i++) // Covariances updates
			{	
				// nur diagonal Elemente berechnen
				for(int d = 0; d < dim; d++)
				{
					sparse_tmsigma[d] += DataSet[i][d]*DataSet[i][d]*pix(i, j);
				}
			}
			
			NICE::Vector sparse_tmsigma2(dim);
			for(int d = 0; d < dim; d++)
			{
				sparse_tmsigma[d] += 1e-5f;
				sparse_tmsigma2[d] = sparse_sigma[j][d]+tmu2[d]*tmu2[d];
				sparse_sigma[j][d] = (sparse_tmsigma[d]+tau*sparse_tmsigma2[d])/(E[j]+tau)-(mu[j][d]*mu[j][d]);
			}				
			sparse_inv_sigma[j] = diagInverse(sparse_sigma[j]);
				
			log_det_sigma[j] = logdiagDeterminant(sparse_sigma[j]);
		}
		if(comp)
		{
			double dist = compare();
			cout << "dist for iteration " << iter << endl;
			cout << "d: " << dist << endl;
		}	
	}
#ifdef DEBUG
	cerr << "GMM::finished EM after " << iter << " iterations" << endl;
#endif
	return iter;
}

int GMM::getBestClass(const NICE::Vector &v, double *bprob)
{
	int bestclass = 0;
	double maxprob = logpdfState(v,0)+log(priors[0]);  // log(P(x|i))
	
	for(int i = 1; i < gaussians; i++)
	{
		double prob = logpdfState(v,i)+log(priors[i]);  // log(P(x|i))
		
		if(prob > maxprob)
		{
			maxprob = prob;
			bestclass = i;
		}
	}
	
	if(bprob != NULL)
		*bprob = maxprob;
	
	return bestclass;
}
		
void GMM::getProbs(const NICE::Vector &vin, SparseVector &outprobs)
{
	outprobs.clear();
	outprobs.setDim(gaussians);
	Vector probs;
	getProbs(vin, probs);
	for(int i = 0; i < gaussians; i++)
	{
		if(probs[i] > 1e-7f )
			outprobs[i] = probs[i];
	}
}
		
void GMM::getProbs(const NICE::Vector &vin, Vector &outprobs)
{
	Vector probs;
	probs.resize(gaussians);
	probs.set(0.0);
		
	double probsum = 0.0;
	double maxp = -numeric_limits<double>::max();
	for(int i = 0; i < gaussians; i++)
	{	
		probs[i] = logpdfState(vin, i) + log(priors[i]);  // log(P(x|i))
		maxp = std::max(maxp, probs[i]);		
	}
		
	for(int i = 0; i < gaussians; i++)
	{
		probs[i] -= maxp;
		probs[i] = exp(probs[i]);
		probsum += probs[i];
	}

	// normalize probs
#pragma omp parallel for
	for(int i = 0; i < gaussians; i++)
	{
		probs[i] /= probsum;
	}
	outprobs = probs;
}

void GMM::getFisher(const NICE::Vector &vin, SparseVector &outprobs)
{
	int size = gaussians*2*dim;
	outprobs.clear();
	outprobs.setDim(size);
	
	int counter = 0;
	
	NICE::Vector classprobs;
	classprobs.resize(gaussians);
	classprobs.set(0.0);

	double maxp = -numeric_limits<double>::max();
	for(int i = 0; i < gaussians; i++)
	{
		classprobs[i] = logpdfState(vin, i) + log(priors[i]);  // log(P(x|i))

		maxp = std::max(maxp, classprobs[i]);
	}

	for(int i = 0; i < gaussians; i++)
	{
		double p = classprobs[i] - maxp;
		
		p = exp(p);
		
		for(int d = 0; d < dim; d++)
		{
			double diff = vin[d]-mu[i][d];
			
			double sigma2 = sparse_sigma[i][d]*sparse_sigma[i][d];
			double sigma3 = sigma2*sparse_sigma[i][d];
			double normmu = sqrt(priors[i]/sigma2);
			double normsig = sqrt((2.0*priors[i])/sigma2);
			
			double gradmu = (p * (diff/sigma2))/normmu;
			double gradsig = (p * ((diff*diff)/sigma3 - 1.0/sparse_sigma[i][d])) /normsig;
			if(fabs(gradmu) > 1e-7f)
				outprobs[counter] = gradmu;
			counter++;
			
			if(fabs(gradsig) > 1e-7f)
				outprobs[counter] = gradsig;
			counter++;
		}
	}
}

void GMM::cluster ( const VVector & features, VVector & prototypes, vector<double> & weights, vector<int> & assignment )
{
	computeMixture(features);

	prototypes.clear();
	weights.clear();
	assignment.clear();

	for(int i = 0; i < (int)features.size(); i++)
	{
		int c = getBestClass(features[i]);
		assignment.push_back(c);
		
		weights.push_back(priors[c]);
	}
	for(int c = 0; c < gaussians; c++)
		prototypes.push_back(mu[c]);
	
	cout << "tau: " << tau << endl;
}

void GMM::saveData(const std::string filename)
{
	ofstream fout( filename.c_str() );

	fout << gaussians << " " << dim << endl;
	
	mu >> fout;
	fout << endl;
			
	for(int n = 0; n < gaussians; n++)
	{
		fout << sparse_sigma[n] << endl;
	}
	
	for(int n = 0; n < gaussians; n++)
	{
		fout << priors[n] << " ";
	}

	fout.close();
}

bool GMM::loadData(const std::string filename)
{
	cerr << "read GMM Data" << endl;
	ifstream fin( filename.c_str() );

	if(fin.fail())
	{
		cerr << "GMM: Datei " << filename << " nicht gefunden!" << endl;
		return false;
	}

	fin >> gaussians;
	fin >> dim;
	cdimval = log(pow(2*3.14159,dim));
	mu.clear();

	for(int i = 0; i < gaussians; i++)
	{
		NICE::Vector tmp;
		fin >> tmp;
		mu.push_back(tmp);
	}
	
	log_det_sigma.clear();

	for(int n = 0; n < gaussians; n++)
	{
		NICE::Matrix _sigma;
		NICE::Vector _sparse_sigma;
		
		_sparse_sigma =NICE::Vector(dim);
		fin >> _sparse_sigma;
		sparse_sigma.push_back(_sparse_sigma);
				
		sparse_inv_sigma.push_back(diagInverse(sparse_sigma[n]));
		log_det_sigma.push_back(logdiagDeterminant(sparse_sigma[n]));
	}

	for(int n = 0; n < gaussians; n++)
	{
		double tmpd;
		fin >> tmpd;
		priors.push_back(tmpd);
	}
	fin.close();

	cerr << "reading GMM Data finished" << endl;
	return true;
}

void GMM::getParams(VVector &mean, VVector &sSigma, vector<double> &p)
{
	mean = mu;
	sSigma.resize(gaussians);
	p.clear();
	for(int i = 0; i < gaussians; i++)
	{
		sSigma[i] = sparse_sigma[i];
		p.push_back(priors[i]);
	}
}

void GMM::setCompareGM(VVector mean, VVector sSigma, vector<double> p)
{
	mu2 = mean;
	sparse_sigma2 = sSigma;
	priors2 = vector<double>(p);
}
		
double GMM::kPPK(NICE::Vector sigma1, NICE::Vector sigma2, NICE::Vector mu1, NICE::Vector mu2, double p)
{
	double d = mu1.size();
	
	double det1 = 1.0;
	double det2 = 1.0;
	double det3 = 1.0;
	double eval = 0.0;
	for(int i = 0; i < d; i++)
	{
		det1 *= sigma1[i];
		det2 *= sigma2[i];
		double sigma = 1.0/((1.0/sigma1[i]+1.0/sigma2[i])*p);
		det3 *= sigma;
		double mu = p*(mu1[i]*sigma1[i]+mu2[i]*sigma2[i]);
		eval += 0.5*mu*sigma*mu - (p*mu1[i]*mu1[i])/(2.0*sigma1[i]) - (p*mu2[i]*mu2[i]) / (2.0*sigma2[i]);
	}
	
	return (pow(2.0*M_PI, ((1.0-2.0*p)*d)/2.0) * pow(det1, -p/2.0) * pow(det2, -p/2.0) * pow(det3, 0.5) * exp(eval));
}

double GMM::compare()
{
	double distkij = 0.0;
	double distkjj = 0.0;
	double distkii = 0.0;
	for(int i = 0; i < gaussians; i++)
	{
		for(int j = 0; j < gaussians; j++)
		{
			double kij = kPPK(sparse_sigma[i],sparse_sigma2[j], mu[i], mu2[j], 0.5);
			double kii = kPPK(sparse_sigma[i],sparse_sigma[j], mu[i], mu[j], 0.5);
			double kjj = kPPK(sparse_sigma2[i],sparse_sigma2[j], mu2[i], mu2[j], 0.5);
			kij = priors[i]*priors2[j]*kij;
			kii = priors[i]*priors[j]*kii;
			kjj = priors2[i]*priors2[j]*kjj;
			distkij += kij;
			distkii += kii;
			distkjj += kjj;
		}
	}
	return (distkij / (sqrt(distkii)*sqrt(distkjj)));
}

void GMM::comparing(bool c)
{
	comp = c;
}