/** 
* @file LocalizationAnalysis.cpp
* @brief Some methods for localization analysis
* @author Erik Rodner
* @date 09/01/2008

*/
#include <iostream>

#include "LocalizationAnalysis.h"

using namespace OBJREC;

using namespace std;
using namespace NICE;

double LocalizationAnalysis::calcAveragePrecision ( const vector<double> & recall,
			      const vector<double> & precision )
{
	// This is the PASCAL VOC version of computing 
	// a kind of interpolated average precision
	// "The intention in interpolating the precision/recall curve
	// in this way is to reduce the impact of the “wiggles” in
	// the precision/recall curve, caused by small variations in the
	// ranking of examples"
	//
	// implementation:
	// \int_0^1 (max_{recall(t') >= t) precision(t')) dt
	//
    double avgp = 0.0;
    for (double t = 0.0 ; t <= 1.0 ; t+=0.1 )
    {
		double p = 0.0;
		for ( uint k = 0 ; k < recall.size() ; k++ )
			if ( recall[k] >= t ) {
				if ( precision[k] >= p ) 
					p = precision[k];
			}

		avgp += p / 11;
    }

    return avgp;
}

double LocalizationAnalysis::calcAveragePrecisionPrecise ( const vector<double> & recall,
			      const vector<double> & precision )
{
	// average precision is defined as:
	// int_0^1 precision(recall) drecall 
	// = int_0^1 precision(t) recall'(t) dt
    double avgp = 0.0;
	double oldrec = recall[0];

	if ( (recall[0] != 0.0) || (recall[recall.size()-1] != 1.0) )
		fthrow(Exception, "Recall/Precision values have to include (0,x) and (1,y)");

    for (uint t = 1 ; t < recall.size(); t++ )
    {
		double prec = precision[t];
		double rec = recall[t];
		// dt of recall' and of the integral are
		// cancelling out
		double drec = rec - oldrec;

		avgp += prec * drec;

		oldrec = rec;
    }

    return avgp;
}

double LocalizationAnalysis::trapezoidArea ( double x1, double y1,
		double x2, double y2 )
{
    return ( x2 - x1 ) * ( y1 + y2 ) / 2.0;
}

double LocalizationAnalysis::calcAreaUnderROC ( const vector<double> & fprate,
		const vector<double> & tprate )
{
	if ( tprate.size() != fprate.size() )
		fthrow(Exception, "True positive rates and false positives rates do not match according to the" <<
			" size of their vectors.");

	double area = 0.0;
	double oldfpr = 0.0;
	double oldtpr = 0.0;

	for ( uint k = 0 ; k < tprate.size() ; k++ )
	{
		double fpr = fprate[k];
		double tpr = tprate[k];
		double newArea = trapezoidArea ( oldfpr, oldtpr, fpr, tpr );	

		area += newArea;
	
		oldfpr = fpr;
		oldtpr = tpr;
	}

	return area;
}

void LocalizationAnalysis::matchResults ( LocalizationResult *l, 
					  const LocalizationResult *l_gt,
					  vector< pair<double, int> > & results,
					  double matchThreshold,
					  bool rejectMultipleDetections )
{
    if ( l == NULL ) return;

    // stat: tp fp tn fn
    set<SingleLocalizationResult *> matchedHypotheses;
    l->sortDescendingConfidence();
    double measureMax = 0.0;

    for ( LocalizationResult::const_iterator i = l->begin();
	    i != l->end(); i++ )
    {
		SingleLocalizationResult *slr_i = *i;
		SingleLocalizationResult *matched_gt = NULL;
		int classno_i = slr_i->r->classno;
		double confidence = slr_i->r->confidence();

		if ( l_gt != NULL )
			for ( LocalizationResult::const_iterator j = l_gt->begin();
			j != l_gt->end(); j++ )
			{
			SingleLocalizationResult *slr_j = *j;

			if ( classno_i == slr_j->r->classno ) {
				double m = slr_i->getBBOverlapMeasure ( *slr_j );
				if ( m > measureMax ) measureMax = m;
				if ( m > matchThreshold )
				{
				if ( matched_gt != NULL )
				{
					fprintf (stderr, "Something strange happened: object hypothesis matched with two different ground truth annotations !! : Overlap issue\n");
					fprintf (stderr, "match score %f\n", m );
					fprintf (stderr, "match threshold %f\n", matchThreshold);
					fprintf (stderr, "sizes algorithm: %zd; sizes groundtruth: %zd\n", l->size(), l_gt->size() );
					int xi, yi, xa, ya;
					slr_i->getBoundingBox ( xi, yi, xa, ya );
					fprintf (stderr, "algorithm: %d %d %d %d\n", xi, yi, xa, ya );
					slr_j->getBoundingBox ( xi, yi, xa, ya );
					fprintf (stderr, "groundtruth: %d %d %d %d\n", xi, yi, xa, ya );
				}
				matched_gt = slr_j;
				}
			}
			}

		if ( matched_gt == NULL )
		{
			results.push_back ( pair<double, int> ( confidence, 0 ) );
		} else {
			if ( matchedHypotheses.find(matched_gt) == matchedHypotheses.end() )
			{
			matchedHypotheses.insert ( matched_gt );
			results.push_back ( pair<double, int> ( confidence, 1 ) );
			} else {
			// Multiple Detections
			if (rejectMultipleDetections)
				results.push_back ( pair<double, int> ( confidence, 0 ) );
			}
		}
    }
}

void LocalizationAnalysis::calcRecallPrecisionCurve ( const vector< pair<double, int> > & results_unsorted,
						      int count_positives,
						      vector<double> & thresholdsv,
						      vector<double> & recallv,
						      vector<double> & precisionv )
{
    vector<pair<double, int> > results;

    for ( vector<pair< double, int > >::const_iterator i = results_unsorted.begin();
	 i != results_unsorted.end(); i++ )
		results.push_back ( pair<double, int> ( - i->first, i->second ) );
    
    sort ( results.begin(), results.end() );
	thresholdsv.clear();
	recallv.clear();
	precisionv.clear();

    long fp = 0;
    long tp = 0;
    for ( vector<pair<double,int> >::const_iterator i = results.begin();
						    i != results.end();
						    i++ )
    {
		double confidence = - i->first;
		int ok = i->second;

		// at this step we have virtually defined a threshold of confidence-epsilon
		tp += ok;
		fp += 1-ok;

		// wikipedia.org
		// Recall in this context is defined as the number of true positives divided by the total number of elements that actually belong 
		// to the positive class (i.e. the sum of true positives and false negatives, which are items which were not labeled as belonging to the positive class but should have been).
		double recall = tp / (double)count_positives;
		// In a statistical classification task, the Precision for a class is the number of true positives 
		// (i.e. the number of items correctly labeled as belonging to the positive class) divided 
		// by the total number of elements labeled as belonging to the positive class (i.e. the sum of true positives and false positives, 
		// which are items incorrectly labeled as belonging to the class). 
		double precision = tp / (double)(fp+tp);

		thresholdsv.push_back ( confidence );
		recallv.push_back ( recall );
		precisionv.push_back ( precision );
    }

	if ( recallv.size() > 0 ) 
	{
		if ( recallv[0] != 0.0 )
		{
			// label everything as negative
			recallv.insert ( recallv.begin(), 0.0 ); 
			precisionv.insert ( precisionv.begin(), 0.0 ); 
			thresholdsv.insert ( thresholdsv.begin(), std::numeric_limits<double>::max() );
		}
		
		if ( recallv[recallv.size() - 1] != 0.0 )
		{
			// label everything as positive
			// fp = count_negatives = results.size() - count_positives
			// tp = count_positives
			// recall = count_positives / count_positives = 1
			// precision = count_positives / ( results.size() - count_positives + count_positives ); 
			//           = count_positives / results.size()
			recallv.insert ( recallv.end(), 1.0 ); 
			precisionv.insert ( precisionv.end(), count_positives / (double)results.size() ); 
			thresholdsv.insert ( thresholdsv.end(), - std::numeric_limits<double>::max() );
		}

	}
		
}

void LocalizationAnalysis::calcROCCurve ( const vector< pair<double, int> > & results_unsorted,
						      int count_positives,
						      int count_negatives,
						      vector<double> & thresholdsv,
						      vector<double> & fprates,
						      vector<double> & tprates )
{
	if ( (count_positives <= 0) )
		fthrow(Exception, "ROC Calculation: number of positive examples is zero");
	if ( (count_negatives <= 0) )
		fthrow(Exception, "ROC Calculation: number of negative examples is zero");

    vector<pair<double, int> > results;

	// sort the results descending
    for ( vector<pair< double, int > >::const_iterator i = results_unsorted.begin();
	 i != results_unsorted.end(); i++ )
		results.push_back ( pair<double, int> ( - i->first, i->second ) );
    
    sort ( results.begin(), results.end() );

    long fp = 0;
    long tp = 0;
	double oldconfidence = - numeric_limits<double>::max();
	// descending order of the results !!
    for ( vector<pair<double,int> >::const_iterator i = results.begin();
						    i != results.end();
						    i++ )
    {
		// everything above the confidence
		// is regarded as positive
		double confidence = - i->first;
		int ok = i->second;

		if ( (ok != 1) && (ok != 0) )
			fthrow(Exception, "ROC Calculation: results should be labeled as 0 or 1 !");

		tp += ok;
		fp += 1-ok;

		//  This test ensures that we wait for the end of a span of
		//  equal probabilities before emitting a point.
		//  Instead of computing a pessimistic version of the ROC curve,
		//  we use the expected ROC curve (see Tom Fawcetts paper)
		//  I changed this in 7/6/2011 to ensure the compability of ROC.pl
		if ( confidence == oldconfidence ) 
			continue;

		double tprate = tp / (double)count_positives;
		double fprate = fp / (double)count_negatives;

		thresholdsv.push_back ( confidence );
		fprates.push_back ( fprate );
		tprates.push_back ( tprate );

		oldconfidence = confidence;
    }

	if ( tprates.size() > 0 ) 
	{
		if ( tprates[0] != 0.0 )
		{
			// label everything as negative
			tprates.insert ( tprates.begin(), 0.0 ); 
			fprates.insert ( fprates.begin(), 0.0 ); 
			thresholdsv.insert ( thresholdsv.begin(), std::numeric_limits<double>::max() );
		}
		
		if ( tprates[tprates.size() - 1] != 1.0 )
		{
			// label everything as positive
			tprates.insert ( tprates.end(), 1.0 ); 
			fprates.insert ( fprates.end(), 1.0 ); 
			thresholdsv.insert ( thresholdsv.end(), - std::numeric_limits<double>::max() );
		}

	}
}

LocalizationAnalysis::LocalizationAnalysis()
{
}

LocalizationAnalysis::~LocalizationAnalysis()
{
}