/** 
* @file FPCBoosting.cpp
* @brief implementation of boosting algorithms 
* @author Erik Rodner
* @date 04/24/2008

*/
#include <iostream>

#include "vislearning/classifier/fpclassifier/randomforest/FPCRandomForests.h"

#include "vislearning/baselib/Gnuplot.h"

#include "FPCBoosting.h"
#include "FPCFullSearch.h"

using namespace OBJREC;

using namespace std;
// refactor-nice.pl: check this substitution
// old: using namespace ice;
using namespace NICE;


#undef DEBUG_BOOST
#define ESTIMATETRAINERROR
#undef DEBUG_ERROR_ESTIMATION

static inline double
pRatio( double val )
{
    const double eps = 1e-5;

    if ( val < eps ) val = eps;
    if ( val > 1 - eps ) val = 1 - eps;
    return val/(1. - val);
}



FPCBoosting::FPCBoosting( const Config *_conf, 
				// refactor-nice.pl: check this substitution
				// old: string section ) : conf(_conf)
				std::string section ) : conf(_conf)
{
    // refactor-nice.pl: check this substitution
    // old: string weakClassifier_s = conf->gS(section, "weak_classifier", "full_search" );
    std::string weakClassifier_s = conf->gS(section, "weak_classifier", "full_search" );

    if ( weakClassifier_s == "random_forest" )
    {
	weakClassifier = new FPCRandomForests ( _conf, "RandomForest" );
	memory_efficient = _conf->gB("RandomForest", "memory_efficient", false );
    } else if ( weakClassifier_s == "full_search" ) {
	fprintf (stderr, "Boost: using full search methods\n");
	weakClassifier = new FPCFullSearch ( _conf );
	memory_efficient = false;
    } else {
	fprintf (stderr, "FPCBoosting: weak classifier type unknown !\n");
	exit(-1);
    }

    // refactor-nice.pl: check this substitution
    // old: string boosting_method_s = conf->gS(section, "method", "realboost" );
    std::string boosting_method_s = conf->gS(section, "method", "realboost" );

    if ( boosting_method_s == "realboost" )
    {
	boosting_method = BOOSTINGMETHOD_REAL_ADABOOST;
    } else if ( boosting_method_s == "adaboost" ) {
	boosting_method = BOOSTINGMETHOD_ADABOOST;
    } else if ( boosting_method_s == "gentleboost" ) {
	boosting_method = BOOSTINGMETHOD_GENTLEBOOST;
    } else {
	fprintf (stderr, "FPCBoosting: boosting method unknown !\n");
	exit(-1);
    }

    classwise_normalization = _conf->gB(section, "classwise_normalization", false );
    positive_class = _conf->gI(section, "positive_class", 1 );
    maxRounds = _conf->gI(section, "max_rounds", 100 );
}

FPCBoosting::~FPCBoosting()
{
    if ( weakClassifier != NULL )
	delete weakClassifier;

    for ( StrongClassifier::iterator i = strongClassifier.begin();
				     i != strongClassifier.end();
				     i++ )
	delete ( i->third );
}


ClassificationResult FPCBoosting::classify ( Example & pce )
{

    FullVector overall_distribution;
    for ( StrongClassifier::const_iterator i = strongClassifier.begin();
					   i != strongClassifier.end();
					   i++ )
    {
	double alpha = i->first;
	double beta  = i->second;
	FeaturePoolClassifier *classifier = i->third;
	ClassificationResult r = classifier->classify ( pce );

	if ( boosting_method == BOOSTINGMETHOD_REAL_ADABOOST )
	{
	    // transform probabilities into nasty scores
	    double p = r.scores[positive_class];
	    r.scores[positive_class] = 0.5 * log ( pRatio ( p ) );
	    r.scores[0] = - r.scores[positive_class];
	}

	if ( overall_distribution.empty() )
	{
	    overall_distribution = r.scores;
	    overall_distribution.multiply(alpha);
	} else {
	    overall_distribution.add ( r.scores, alpha );
	}
	overall_distribution.add ( beta );
    }

    int overall_classno = 0;
    overall_classno = overall_distribution.maxElement();

    return ClassificationResult ( overall_classno, overall_distribution );
}
	
void FPCBoosting::train ( FeaturePool & fp,
      Examples & examples )
{
    maxClassNo = examples.getMaxClassNo();

    FeatureStorage featureStorage;

#ifdef DEBUG_BOOST
    fprintf (stderr, "FPCBoosting : training examples %d\n", (int)examples.size() );
#endif

    boosting ( featureStorage, fp, examples );
}

void FPCBoosting::normalizeWeights ( Examples & examples ) const
{
    if ( classwise_normalization )
    {
	double sump = 0.0;
	double sumn = 0.0;
	int np = 0;
	int nn = 0;

	for ( Examples::const_iterator i = examples.begin();
					       i != examples.end();
					       i++ )
	    if ( i->first == positive_class )
	    {
		sump += i->second.weight;
		np++;
	    } else {
		sumn += i->second.weight;
		nn++;
	    }

	if ( fabs(sump) < 1e-10 ) sump = np;
	if ( fabs(sumn) < 1e-10 ) sump = nn;
	
	for ( Examples::iterator i = examples.begin();
					       i != examples.end();
					       i++ )
	    if ( i->first == positive_class )
		i->second.weight /= 2*sump;
	    else
		i->second.weight /= 2*sumn;


    } else {
	double sum = 0.0;
	for ( Examples::const_iterator i = examples.begin();
					       i != examples.end();
					       i++ )
	    sum += i->second.weight;

	if ( fabs(sum) < 1e-10 ) sum = examples.size();
	for ( Examples::iterator i = examples.begin();
					       i != examples.end();
					       i++ )
	    i->second.weight /= sum;
    }

}

void FPCBoosting::boosting ( const FeatureStorage & featureStorage,
				FeaturePool & fp,
				Examples & examples )
{
    normalizeWeights ( examples );
#ifdef ESTIMATETRAINERROR
    vector<double> error_boosting;
    vector<double> weak_error_boosting;
#endif

    for ( uint iteration = 0 ; iteration < (uint)maxRounds ; iteration++ )
    {
	// fit a classifier using the current weights
	FeaturePoolClassifier *best = weakClassifier->clone();
	best->train ( fp, examples );

	// ---------------------------------------------
	// Estimate the error of the current hypothesis
	// ---------------------------------------------
	double error = 0.0;
	long int index = 0;
	vector<double> h_values (examples.size(), 0.0);
	vector<int> h_results (examples.size(), 0);
	
	for ( Examples::iterator i = examples.begin();
					   i != examples.end();
					   i++, index++ )
	{
	    double weight = i->second.weight;

	    ClassificationResult r = best->classify ( i->second );
	    double p_value = r.scores[positive_class];
	    int h_result = r.classno == positive_class ? 1 : -1;
	    double h_value;
	    
	    // assume p_value is between [0:1]
	    if ( boosting_method == BOOSTINGMETHOD_REAL_ADABOOST )
	    {
		if ( (p_value < 0.0) || (p_value > 1.0) )
		{
		    fprintf (stderr, "FPCBoosting: do not use real adaboost with hypothesis values outside of [0,1]\n");
		    exit(-1);
		}
		h_value = 0.5 * log (pRatio ( p_value ));
	    } else {
		h_value = p_value;
	    }

	    assert ( index < (int)h_values.size() );
	    assert ( index < (int)h_results.size() );
	    h_values[index] = h_value;
	    h_results[index] = h_result;

	    if ( r.classno != i->first )
		error += weight;
	   
	    #ifdef DEBUG_BOOST
	    fprintf (stderr, "FPCBoosting:w(%ld) = %f gt %d est %d\n", index, weight, i->first, h_result );
	    #endif

	    if ( memory_efficient ) i->second.ce->dropPreCached();
	}

	#ifdef DEBUG_ERROR_ESTIMATION
	fprintf (stderr, "Boost: iteration %zd error %lf\n", iteration, error );
	FPCFullSearch *search = dynamic_cast< FPCFullSearch * > ( best );
	
	/*
	if ( fabs(error - search->last_error) > 1e-5 ) {
	    fprintf (stderr, "Boost: FIX THIS BUG postest=%lf preest=%lf %e\n", error, search->last_error, fabs(error - search->last_error) );
	    exit(-1);
	}
	*/
	
	if ( error > 0.5 )
	{
	    fprintf (stderr, "Boost: weak hypothesis with an error > 0.5 ? postest=%lf preeest=%lf\n", error, search->last_error);
	    exit(-1);
	}
	#endif

	double likelihood_ratio = 1.0;
	if (boosting_method == BOOSTINGMETHOD_REAL_ADABOOST) 
	{
	    strongClassifier.push_back ( triplet<double, double, FeaturePoolClassifier *> ( 1.0, 0.0, best ) );
	} else if ( (boosting_method == BOOSTINGMETHOD_GENTLEBOOST) ) {
	    strongClassifier.push_back ( triplet<double, double, FeaturePoolClassifier *> ( 1.0, 0.0, best ) );
	} else {
	    // likelihood_ratio corresponds to \beta_t in the 
	    // Viola&Jones Paper
	    likelihood_ratio = pRatio ( error );
	    double alpha = - log ( likelihood_ratio );
	    double beta  = alpha/2;

	    #ifdef DEBUG_BOOST
	    fprintf (stderr, "estimated parameters: lratio=%f alpha=%f beta=%f\n", likelihood_ratio, alpha, beta);
	    #endif
	    strongClassifier.push_back ( triplet<double, double, FeaturePoolClassifier *> ( alpha, beta, best ) );
	}

	#ifdef ESTIMATETRAINERROR
	    // --------- estimate training error
	    double error_sum = 0.0;
	    index = 0;
	    for ( Examples::iterator i = examples.begin();
					       i != examples.end();
					       i++, index++ )
	    {
		ClassificationResult r = classify ( i->second );

		if ( r.classno != i->first )
		    error_sum += 1.0 / examples.size();
	    }
	    fprintf (stderr, "Boost: training error %f\n", error_sum );
	    
	    error_boosting.push_back ( error_sum );
	    weak_error_boosting.push_back ( error );
	#endif

	// ---------- readjust weights
	index = 0;
	for ( Examples::iterator i = examples.begin();
					   i != examples.end();
					   i++, index++ )
	{
	    double weight = i->second.weight;
	    assert ( index < (int)h_values.size() );
	    assert ( index < (int)h_results.size() );
	    double h_value = h_values[index];
	    int h_result = h_results[index];
	    int y = (i->first == positive_class) ? 1 : -1;
    
	    if (boosting_method == BOOSTINGMETHOD_REAL_ADABOOST)
	    {
		// weight update of Real-AdaBoost
		// as presented by Friedman et al.
		weight *= exp( - y*h_value );
	    } else if (boosting_method == BOOSTINGMETHOD_GENTLEBOOST) {
		// h_value is still between [0,1]
		weight *= exp( - y*h_value );
	    } else {
		// standard AdaBoost weight update
		// according to Viola & Jones
		if ( y == h_result )
		    weight *= likelihood_ratio;
	    }


	    #ifdef DEBUG_BOOST
	    fprintf (stderr, "Boost: iteration %d y %d p %f w %f\n", iteration, y, h_value, weight );
	    #endif

	    i->second.weight = weight;
	}
	normalizeWeights ( examples );
    }

#ifdef ESTIMATETRAINERROR
    Gnuplot gp ("lines");
    vector<double> upper_bound;

    // Compute the upper bound on the training error
    // the formula and an explanation can be found in the
    // paper of Freund & Shapire
    // 2^T can be incorporated into the product
    // and a quick glance at the formula shows that the
    // upper bound is at least below 1 (otherwise this
    // bound we be cound of useless
    double prod = 1.0;
    for ( vector<double>::const_iterator i = weak_error_boosting.begin();
	    i != weak_error_boosting.end(); i++ )
    {
	double epsilon_t = *i;
	prod *= 2 * sqrt( epsilon_t * (1 - epsilon_t) );
	// according to the improvment of the upper bound in the paper
	upper_bound.push_back ( 0.5*prod );
    }

    gp.plot_x ( error_boosting, "Boosting Training Error" );
    gp.plot_x ( weak_error_boosting, "Error of the weak learner" );
    gp.plot_x ( upper_bound, "Upper bound of the training error (AdaBoost using a Soft-Decision)" );

    #ifndef NOVISUAL
    // refactor-nice.pl: check this substitution
    // old: GetChar();
    getchar();
    #else
    getchar();
    #endif
#endif
}


void FPCBoosting::restore (istream & is, int format)
{
    // refactor-nice.pl: check this substitution
    // old: string tag;
    std::string tag;

    weakClassifier->maxClassNo = maxClassNo;

    while ( (is >> tag) && (tag == "WEAKCLASSIFIER") )
    {
	FeaturePoolClassifier *classifier = weakClassifier->clone();

	double alpha;
	double beta;

	is >> alpha;
	is >> beta;
	classifier->restore ( is );

	strongClassifier.push_back ( triplet<double, double, FeaturePoolClassifier *>
	    ( alpha, beta, classifier ) );
    }
    cerr << "TAG: " << tag << endl;

    assert ( strongClassifier.size() > 0 );

}

void FPCBoosting::store (ostream & os, int format) const
{
    for ( StrongClassifier::const_iterator i = strongClassifier.begin();
				           i != strongClassifier.end();
					   i++ )
    {
	const FeaturePoolClassifier *classifier = i->third;
	double alpha = i->first;
	double beta  = i->second;
    
	os << "WEAKCLASSIFIER" << endl;

	os << alpha << endl;
	os << beta << endl;

	classifier->store (os);

	os << "ENDWEAKCLASSIFIER" << endl;

    }
}

void FPCBoosting::clear ()
{
    strongClassifier.clear();
}

FeaturePoolClassifier *FPCBoosting::clone () const
{
    FPCBoosting *fpcBoost = new FPCBoosting ();

    fpcBoost->maxRounds = maxRounds;
    fpcBoost->weakClassifier = weakClassifier->clone();
    fpcBoost->positive_class = positive_class;
    fpcBoost->memory_efficient = memory_efficient;
    fpcBoost->maxClassNo = maxClassNo;

    return fpcBoost;

}

void FPCBoosting::setComplexity ( int size )
{
    fprintf (stderr, "FPCBoosting: set complexity to %d, overwriting current value %d\n", 
	size, maxRounds );
    maxRounds = size;
}