#ifdef NOVISUAL
#warning "testKMeans needs ICE with visualization !!"
    int main (int argc, char **argv) {};
#else

#ifndef M_PI
#define M_PI 3.14159265358979323846f
#endif

#include <fstream>
#include "core/vector/VectorT.h"
#include "core/vector/MatrixT.h"
#include "core/image/ImageT.h"
#include "core/imagedisplay/ImageDisplay.h"
#include <iostream>
#include <core/image/CrossT.h>
#include "vislearning/math/cluster/GMM.h"


using namespace std;
using namespace NICE;
using namespace OBJREC;
		 
int main (int argc, char **argv)
{
	Config conf ( argc, argv );
	
	std::string params = conf.gS("main", "params");
	std::string featfile = conf.gS("main", "feats");
	
	ifstream fin( params.c_str() );
	
	if(fin.bad())
	{
		cout << "there were errors while opening the file " << params << endl;
	}
	
	int dim, gaussians, feats;
	double tmp;
	fin >> dim;
	for(int i = 1; i < dim; i++)
	{
		fin >> tmp;
	}
	fin >> gaussians;
	for(int i = 1; i < dim; i++)
	{
		fin >> tmp;
	}
	fin >> feats;
	for(int i = 1; i < dim; i++)
	{
		fin >> tmp;
	}
	
	vector<double> prior1;
	
	VVector mean1(gaussians, dim);
	VVector sigma1(gaussians, dim);
	
	for(int i = 0; i < gaussians; i++)
	{
		fin >> tmp;
		prior1.push_back(tmp);
		for(int d = 1; d < dim; d++)
		{
			fin >> tmp;
		}
		for(int d = 0; d < dim; d++)
		{
			fin >> mean1[i][d];
		}
		for(int d1 = 0; d1 < dim; d1++)
		{
			for(int d2 = 0; d2 < dim; d2++)
			{
				if(d1 == d2)
					fin >> sigma1[i][d1];
				else
					fin >> tmp;
			}
		}
	}
	
	int featsize;
	fin >> featsize;
	fin.close();
	
	ifstream fin2( featfile.c_str() );
	if(fin2.bad())
	{
		cout << "there were errors while opening the file " << featfile << endl;
	}
	
	VVector x;
	VVector prototypes;
	std::vector<double> weights;
	std::vector<int> assignments;
	
	for(int i = 0; i < featsize; i++)
	{
		NICE::Vector f(dim);
		for(int d = 0; d < dim; d++)
		{
			fin2 >> f[d];
		}
		x.push_back(f);
	}

	GMM clusteralg( &conf, gaussians);
#if 1
	//clusteralg->computeMixture ( x );
	clusteralg.setGMMtoCompareWith(mean1, sigma1, prior1);
	clusteralg.setCompareTo2ndGMM(true);
	clusteralg.cluster ( x, prototypes, weights, assignments );
	
	double dist = clusteralg.compareTo2ndGMM();
	
	cout << "dist1: " << endl << dist << endl;
	
	int width = 500;
	int height = 500;
	NICE::Image panel (width, height);
	NICE::Image overlay (width, height);
	panel.set(255);
	overlay.set(0);
	width--;
	height--;
	for ( size_t i = 0 ; i < prototypes.size() ; i++ )
	{
	    // refactor-nice.pl: check this substitution
		fprintf (stderr, "prototype %d (%f, %f)\n", (int)i, prototypes[i][0], prototypes[i][1] );
		Cross cross ( Coord( (int)(prototypes[i][0]*(double)width), (int)(prototypes[i][1]*(double)height) ), 3 );
		overlay.draw ( cross, i+1 );
	}
	fprintf (stderr, "assignments: %ld\n", assignments.size() );

	for ( size_t i = 0 ; i < assignments.size() ; i++ )
	{
		Cross cross ( Coord( (int)(x[i][0]*(double)width), (int)(x[i][1]*(double)height)) , 3 );
		overlay.draw ( cross, assignments[i]+1 );

	}
	NICE::showImageOverlay ( panel, overlay, "Clustering results" );
	
#else
	double scale = 2.0;
	vector<double> prior2 = prior1;
	VVector sigma2 = sigma1;
	VVector mean2 = mean1;
	for(int i = 0; i < sigma2.size(); i++)
	{
		for(int j = 0; j < sigma2[i].size(); j++)
		{
			sigma2[i][j] /= scale;
		}
	}
	for(int i = 0; i < mean2.size(); i++)
	{
		for(int j = 0; j < mean2[i].size(); j++)
		{
			mean2[i][j] /= scale;
		}
	}
	for(int i = 0; i < prior2.size(); i++)
	{
		prior2[i] /= scale;
	}
	double distkii = 0.0;
	double distkjj = 0.0;
	double distkij = 0.0;
	double dist = 0.0;
	for(int i = 0; i < gaussians; i++)
	{
		for(int j = 0; j < gaussians; j++)
		{
			double kij = clusteralg.kPPK(sigma1[i],sigma2[j], mean1[i], mean2[j], 0.5);
			double kii = clusteralg.kPPK(sigma1[i],sigma1[j], mean1[i], mean1[j], 0.5);
			double kjj = clusteralg.kPPK(sigma2[i],sigma2[j], mean2[i], mean2[j], 0.5);
			kij = prior1[i]*prior2[j]*kij;
			kii = prior1[i]*prior1[j]*kii;
			kjj = prior2[i]*prior2[j]*kjj;
			distkii += kii;
			distkjj += kjj;
			distkij += kij;
		}
	}
	
	double dist2 = distkij / (sqrt(distkii)*sqrt(distkjj));
	cout << "dist: " << dist << endl;
	cout << "dist2: " << dist2 << endl;
	for(int k = 0; k < 10; k++)
	{
		dist = 0.0;
		for(int i = 0; i < gaussians; i++)
		{
			int j = i+k;
			if(j >= gaussians)
				j -= gaussians;
			double kij = clusteralg.kPPK(sigma1[i],sigma2[j], mean1[i], mean2[j], 0.5);
			kij = prior1[i]*prior2[j]*kij;
			double kii = clusteralg.kPPK(sigma1[i],sigma1[i], mean1[i], mean1[i], 0.5);
			kii = prior1[i]*prior1[j]*kii;
			double kjj = clusteralg.kPPK(sigma2[j],sigma2[j], mean2[j], mean2[j], 0.5);
			kjj = prior2[i]*prior2[j]*kjj;
			double val = kij / (sqrt(kii)*sqrt(kjj));
			if(kii == 0.0 || kjj == 0.0)
				continue;
			dist+=val;
			//printf ("%4.6f ", val*10.0);
			//cout << endl;
		}
		cout << "dist" << 3+k << " " <<  dist << endl;
	}
#endif
    return 0;
}

#endif