/*
    Copyright (C) 2012 Bartosz Broda, Jacek Bochenek, Łukasz Biłeńkij, 
    Damian Jaworski, Wojciech Mazur, Dominik Piasecki, Maciej Piasecki, 
    Radosław Ramocki
    
    This program is free software; you can redistribute it and/or modify it
under the terms of the GNU Lesser General Public License as published by the Free
Software Foundation; either version 3 of the License, or (at your option)
any later version.

    This program is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE. 

    See the LICENSE and COPYING files for more details.
*/

#include "Evaluator.h"

namespace evaluator
{
	using namespace std;
	using namespace smartcomparator::compare;

	Evaluator::Evaluator()
	{

	}

	void Evaluator::InitEvaluator(const std::string &clutoFile, const std::string &classLabelsFile, bool silent)
	{
		_rows = _clusteredRows = 0;
		this->silent = silent;
		_classes.clear();
		_classesWOutliers.clear();
		_classLabels.clear();
		_clusters.clear();
		_clusterIDs.clear();
		ReadClutoFile(clutoFile);
		if (classLabelsFile != "")
			ReadClassLabelsFile(classLabelsFile);
	}

	void Evaluator::ReadClutoFile(const std::string &clutoFile)
	{
		ifstream in_stream(clutoFile.c_str());
		if(!in_stream)
			throw runtime_error("Error reading from file " + clutoFile);

		string line;
		size_t row = 0;

		if (!silent)
			cout << "Reading cluto file" << endl;
		while (getline(in_stream, line, '\n'))
		{
			try
			{
				int space = line.find(" ");
				if (space > -1)
					line = line.substr(0, space);
				int cluster = boost::lexical_cast<int>(line);

				_clusterIDs.push_back(cluster);
				std::map<int, std::vector<size_t> >::iterator iter = _clusters.find(cluster);

				if (iter == _clusters.end()) //initialize vector
					_clusters.insert(make_pair(cluster, vector<size_t>(1, row)));
				else
					iter->second.push_back(row);
			}
			catch (boost::bad_lexical_cast)
			{
				throw runtime_error("Error in cluto file, line number: " + (row+1));
			}
			row++;
		}
		_rows = row;
		if (!silent)
			cout << "Done" << endl;
	}

	void Evaluator::ReadClassLabelsFile(const std::string &classLabelsFile)
	{
		ifstream in_stream(classLabelsFile.c_str());
		if(!in_stream)
			throw runtime_error("Error reading from file " + classLabelsFile);

		string line;
		size_t row = 0;

		if (!silent)
			cout << "Reading class labels file" << endl;
		while (getline(in_stream, line, '\n'))
		{
			_classLabels.push_back(line);

			//put into class->classMembers map (even if it's an outlier)
			std::map<std::string, std::vector<size_t> >::iterator iter = _classesWOutliers.find(line);
			if (iter == _classesWOutliers.end()) //initialize vector
			    _classesWOutliers.insert(make_pair(line, vector<size_t>(1, row)));
			else
			    iter->second.push_back(row);

			//put into class->classMembers map (if it isn't an outlier)
			if (_clusterIDs[row] >= 0) //check
			{
			    std::map<std::string, std::vector<size_t> >::iterator iter2 = _classes.find(line);
			    if (iter2 == _classes.end()) //initialize vector
				_classes.insert(make_pair(line, vector<size_t>(1, row)));
			    else
				iter2->second.push_back(row);
			}
			row++;
		}
		if (_classLabels.size() != _rows)
			throw runtime_error("Class labels file has different number of rows than cluto file - classlabels: " + boost::lexical_cast<string>(_classLabels.size()) + " vs rows: " + boost::lexical_cast<string>(_rows));
		if (!silent)
			cout << "Done" << endl;
	}

	bool value_comparer(map<string, size_t>::value_type &i1, map<string, size_t>::value_type &i2)
	{
		return i1.second < i2.second;
	}

	double Evaluator::CountPurity()
	{
		if (this->_classes.empty())
			throw runtime_error("Class labels not available");
		if (!silent)
			cout << "Counting purity measure" << endl;
		double result = 0.0;
		map<int, vector<size_t> >::iterator cluster;
		_clusteredRows = _rows;

		for (cluster = _clusters.begin(); cluster != _clusters.end(); ++cluster)
		{
			if (cluster->first < 0)
			{
				_clusteredRows -= cluster->second.size();
				if (!silent)
					cout << "Unclustered: " << cluster->second.size() << endl;
				continue;
			}
			vector<size_t> & elems = cluster->second;
			map<string, size_t> labels;
			for (vector<size_t>::iterator elem = elems.begin(); elem != elems.end(); ++elem)
			{
				string label = _classLabels[*elem];
				if (labels.find(label) == labels.end())
					labels.insert(make_pair(label, 1));
				else
					labels[label]++;
			}
			map<string, size_t>::iterator max = max_element(labels.begin(), labels.end(), value_comparer);
			if (!silent)
				cout << "Cluster: " << cluster->first << " label: " << max->first << " elems: " << max->second << endl;
			result += max->second;
		}
		if (!silent)
		{
			cout << "Clustered rows: " << _clusteredRows << endl;
			cout << "Done" << endl;
		}
		return (result / _clusteredRows);
	}

	double Evaluator::CountNMI()
	{
		if (this->_classes.empty())
			throw runtime_error("Class labels not available");
		if (!silent)
			cout << "Counting nmi measure" << endl;
		double result = 0.0;
		map<int, vector<size_t> >::iterator clusterIter;
		map<string, vector<size_t> >::iterator classIter;
		_clusteredRows = _rows;

		//count mutual information - I(Om, C), see equation 16.4
		for (clusterIter = _clusters.begin(); clusterIter != _clusters.end(); ++clusterIter)
		{
			if (clusterIter->first < 0)
			{
				_clusteredRows -= clusterIter->second.size();
				continue;
			}
			for (classIter = _classes.begin(); classIter != _classes.end(); ++classIter)
			{
				size_t commonElems = 0; //number of elems belonging to current cluster & current class at once
				vector<size_t> & elems = clusterIter->second;
				for (vector<size_t>::iterator elem = elems.begin(); elem != elems.end(); ++elem)
				{
					if (_classLabels[*elem] == classIter->first)
						commonElems++;
				}
				if (commonElems == 0) //check it
					continue;
				double logNo = ((double)(_clusteredRows * commonElems)) / ((double)(clusterIter->second.size() * classIter->second.size()));
				result += ((double)commonElems / (double)_clusteredRows) * log(logNo);
			}
		}
		if (!silent)
			cout << "I(Om, C): " << result << endl;

		//count clusters entropy - H(Om), see equation 16.6
		double clustersEntropy = 0.0;
		for (clusterIter = _clusters.begin(); clusterIter != _clusters.end(); ++clusterIter)
		{
			if (clusterIter->first < 0)
				continue;
			//cout << "Elems in cluster: " << clusterIter->first << " " << clusterIter->second.size() << "/" << clusteredRows << endl;
			double number = (double)clusterIter->second.size() / (double)_clusteredRows;
			//cout << "Log number cluster: " << number << endl;
			clustersEntropy += number * log(number);
		}
		clustersEntropy *= -1;
		if (!silent)
			cout << "H(Om):" << clustersEntropy << endl;

		//count classes entropy - H(C), no equation as an example
		double classEntropy = 0.0;
		for (classIter = _classes.begin(); classIter != _classes.end(); ++classIter)
		{
			//cout << "Elems in class: " << classIter->first << " " << classIter->second.size() << "/" << clusteredRows << endl;
			double number = (double)classIter->second.size() / (double)_clusteredRows;
			//cout << "Log number class: " << number << endl;
			classEntropy += number * log(number);
		}
		classEntropy *= -1;
		if (!silent)
			cout << "H(C): " << classEntropy << endl << "Clustered rows: " << _clusteredRows << endl << "Done" << endl;

		return result / ((classEntropy + clustersEntropy)/2);
	}

	size_t CountNewton2(size_t param)
	{
		return param * (param-1) / 2;
	}

	Stats Evaluator::CountStats()
	{
		Stats s;
		s.TP = s.FP = s.TN = s.FN = 0;
		size_t sumTPFP = 0;

		map<int, vector<size_t> >::iterator clusterIter;
		for (clusterIter = _clusters.begin(); clusterIter != _clusters.end(); ++clusterIter)
		{
			if (clusterIter->first < 0)
				continue;
			sumTPFP += CountNewton2(clusterIter->second.size());

			map<string, size_t> clusterElems; //map: (class name, number of class elems in current cluster)
			vector<size_t>::iterator elem = clusterIter->second.begin();
			for (; elem != clusterIter->second.end(); ++elem)
			{
				string currentLabel = _classLabels[*elem];
				if (clusterElems.find(currentLabel) == clusterElems.end())
					clusterElems.insert(make_pair(currentLabel, 1));
				else
					clusterElems[currentLabel]++;
			}
			for (map<string, size_t>::iterator iter = clusterElems.begin(); iter != clusterElems.end(); ++iter)
			{
				s.TP += CountNewton2(iter->second);
			}
		}
		s.FP = sumTPFP - s.TP;

		map<string, vector<size_t> >::iterator classIter;
		for (classIter = _classes.begin(); classIter != _classes.end(); ++classIter)
		{
			s.FN += CountNewton2(classIter->second.size());
		}
		s.FN -= s.TP;

		size_t all = CountNewton2(_rows);

		s.TN = all - s.FN - s.FP - s.TP;

		return s;
	}

	double Evaluator::CountRandIndex()
	{
		if (this->_classes.empty())
			throw runtime_error("Class labels not available");
		if (!silent)
			cout << "Counting Rand Index measure" << endl;
		Stats stats = CountStats();

		if (!silent)
		{
			cout << "TP: " << stats.TP << endl;
			cout << "FP: " << stats.FP << endl;
			cout << "TN: " << stats.TN << endl;
			cout << "FN: " << stats.FN << endl;

			cout << "Done" << endl;
		}
		return (double)(stats.TP + stats.TN) / (double)(stats.TP + stats.TN + stats.FP + stats.FN);
	}

	double Evaluator::CountPrecision(Stats &stat)
	{
		return (double)stat.TP / (double)(stat.TP + stat.FP);
	}

	double Evaluator::CountRecall(Stats &stat)
	{
		return (double)stat.TP / (double)(stat.TP + stat.FN);
	}

	double Evaluator::CountFmeasure(double beta)
	{
		if (this->_classes.empty())
			throw runtime_error("Class labels not available");
		if (!silent)
			cout << "Counting Fmeasure(" << beta << ")" << endl;
		Stats stats = CountStats();

		double precision = CountPrecision(stats);
		double recall = CountRecall(stats);
		double square = beta * beta;
		if (!silent)
		{
			cout << "Precision: " << precision << endl;
			cout << "Recall: " << recall << endl;
			cout << "Done" <<  endl;
		}
		return ((square + 1) * precision * recall) / (square * precision + recall);
	}

	size_t Evaluator::GetRows() const
	{
		return _rows;
	}

	size_t Evaluator::GetClusteredRows() const
	{
		return _clusteredRows;
	}

	size_t Evaluator::GetFoundClusters() const
	{
		size_t size = _clusters.size();
		if (_clusters.find(-1) != _clusters.end())
			size--;
		return size;
	}

	double Evaluator::CountBCubed(double alpha, double & precision, double & recall)
	{
		if (this->_classes.empty())
			throw runtime_error("Class labels not available");
		double precisionSum = 0;
		double recallSum = 0;

		for (size_t i = 0; i < _clusterIDs.size(); ++i)
		{
			int currentClusterId = _clusterIDs[i];
			if (currentClusterId < 0) //omit the outliers
				continue;
			std::vector<size_t> & currentCluster = _clusters[currentClusterId];
			size_t correctnessSum = 0;

			for (size_t j = 0; j < currentCluster.size(); ++j) //count elements from currentCluster with the same class_label
			{
				if (_classLabels[i] == _classLabels[currentCluster[j]])
					correctnessSum++;
			}

			precisionSum += (double)correctnessSum / (double)currentCluster.size();
			std::string currentClassLabel = _classLabels[i];
			recallSum += (double)correctnessSum / (double)_classesWOutliers[currentClassLabel].size();
		}

		//count the averages
		precision = precisionSum / _clusterIDs.size();
		recall = recallSum / _clusterIDs.size();

		if (!silent)
		{
			cout << "BCubed precision: " << precision << endl;
			cout << "BCubed recall: " << recall << endl;
		}

		double F_rev = alpha * (1 / precision) + (1 - alpha) * (1 / recall);

		return (1 / F_rev);
	}

	//---------------------------------------------------------------------------------------------------------------
	//----------------------------------Internal Evaluation Methods--------------------------------------------------
	//---------------------------------------------------------------------------------------------------------------

	double Evaluator::CountSimilarity(size_t i, size_t j)
	{
		double sim;

		if (_sim_matrix)
		{
			sim = (*_sim_matrix->GetRow(i))[j];
		}
		else
		{
			SimilarityMap::iterator iter = _sim_map.find(make_pair(i, j));
			if(iter == _sim_map.end())
			{
				iter = _sim_map.find(make_pair(j, i));
				if(iter == _sim_map.end())
				{
					sim = _cmp_fun->Compare(i, j);
					_sim_map.insert(std::make_pair( std::make_pair(i, j), sim ));
				}
				else
					sim = iter->second;
			}
			else
				sim = iter->second;
		}
		return sim;
	}

	void Evaluator::SetCompareFunction(smartcomparator::compare::CompareFunctionPtr cmp_fun)
	{
		_cmp_fun = cmp_fun;
	}

	void Evaluator::SetSimilarityMatrix(Matrices::SuperMatrixPtr sim_matrix)
	{
		_sim_matrix = sim_matrix;
	}

	double Evaluator::CountDunnIndex(smartcomparator::compare::CompareFunctionPtr cmp_fun, bool useImprovedFunctions, double maxValue)
	{
		double minGlobalDist = maxValue;
		double maxGlobalDiameter = cmp_fun->MinValue();

		std::map<int, std::vector<size_t> >::iterator iter1 = _clusters.begin();

		for (; iter1 != _clusters.end(); ++iter1)
		{
			std::map<int, std::vector<size_t> >::iterator iter2(iter1);
			iter2++;
			for (; iter2 != _clusters.end(); ++iter2)
			{
				double minDist;
				if (useImprovedFunctions)
					minDist = DistanceBetweenClusters2(iter1->second, iter2->second);
				else
					minDist = DistanceBetweenClusters(iter1->second, iter2->second, maxValue);
				if (minDist < minGlobalDist)
					minGlobalDist = minDist;
			}
			double maxDiameter;
			if (useImprovedFunctions)
				maxDiameter = ClusterDiameter2(iter1->second);
			else
				maxDiameter = ClusterDiameter(iter1->second);
			if (maxDiameter > maxGlobalDiameter)
				maxGlobalDiameter = maxDiameter;
		}
		return (minGlobalDist / maxGlobalDiameter);
	}

	double Evaluator::DistanceBetweenClusters(const std::vector<size_t> & vec1, const std::vector<size_t> & vec2, double maxValue)
	{
		double minDist = maxValue;
		for (size_t i = 0; i < vec1.size(); ++i)
		{
			for (size_t j = 0; j < vec2.size(); ++j)
			{
				double dist = 1.0 - CountSimilarity(vec1[i], vec2[j]);
				if (dist < minDist)
					minDist = dist;
			}
		}
		return minDist;
	}

	double Evaluator::ClusterDiameter(const std::vector<size_t> & vec)
	{
		double maxDiameter = _cmp_fun->MinValue();
		for (size_t i = 0; i < vec.size(); ++i)
		{
			for (size_t j = i + 1; j < vec.size(); ++j)
			{
				double dist = 1.0 - CountSimilarity(vec[i], vec[j]);
				if (dist > maxDiameter)
					maxDiameter = dist;
			}
		}
		return maxDiameter;
	}

	double Evaluator::DistanceBetweenClusters2(const std::vector<size_t> &vec1, const std::vector<size_t> &vec2)
	{
		double sumDist = 0.0;
		for (size_t i = 0; i < vec1.size(); ++i)
		{
			for (size_t j = 0; j < vec2.size(); ++j)
			{
				double dist = 1.0 - CountSimilarity(vec1[i], vec2[j]);
				sumDist += dist;
			}
		}
		double sizes = vec1.size() * vec2.size();
		return (1/sizes) * sumDist;
	}

	double Evaluator::ClusterDiameter2(const std::vector<size_t> & vec)
	{
		using namespace Matrices;

		CompressedVectorPtr meanVector = CompressedVectorPtr(new CompressedVector(_cmp_fun->GetMatrix()->NumberOfColumns()));

		for (size_t i = 0; i < vec.size(); ++i)
		{
			CompressedVectorPtr currentRow = _cmp_fun->GetMatrix()->GetRow(vec[i]);
			*meanVector += *currentRow;
		}
		*meanVector /= vec.size();

		double sum = 0.0;
		for (size_t i = 0; i < vec.size(); ++i)
		{
			double dist = 1 - _cmp_fun->ExternalCompare(vec[i], meanVector);
			sum += dist;
		}
		return 2 * (sum / vec.size());
	}

	double Evaluator::CountExpectedDensity()
	{
		double vertices = _clusterIDs.size(); // |V|
		double wG = vertices; // w(G) = |V| + sum_of_similarities
		for (size_t i = 0; i < vertices; ++i)
			for (size_t j = i + 1; j < vertices; ++j)
				wG += CountSimilarity(i, j);

		double theta = log(wG)/log(vertices); // \Theta

		std::map<int, std::vector<size_t> >::iterator cluster = _clusters.begin();

		double expected_density = 0.0;
		for (; cluster != _clusters.end(); ++cluster)
		{
			double vertices_i = cluster->second.size(); // |V_i|
			double wG_i = vertices_i; // w(G_i)
			for (size_t i = 0; i < vertices_i; ++i)
				for (size_t j = i + 1; j < vertices_i; ++j)
					wG_i += CountSimilarity(cluster->second[i], cluster->second[j]);

			double nominator = vertices_i * wG_i;
			double denominator = vertices * pow(vertices_i, theta);
			expected_density += nominator / denominator;
		}
		return expected_density;
	}
}
