/** 
* @file RegGaussianProcess.cpp
* @brief Regression using Gaussian Processes
* @author Erik Rodner
* @date 12/09/2009

*/
#include <iostream>
#include <sstream>

#include "core/optimization/gradientBased/OptimizationAlgorithmFirst.h"
#include "core/optimization/gradientBased/FirstOrderTrustRegion.h"
#include "core/optimization/gradientBased/FirstOrderRasmussen.h"

#include "core/vector/Algorithms.h"
#include "core/algebra/CholeskyRobust.h"
#include "core/algebra/CholeskyRobustAuto.h"

#include "RegGaussianProcess.h"
#include "GPRegressionOptimizationProblem.h"

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



RegGaussianProcess::RegGaussianProcess( const Config *conf, Kernel *kernelFunction, const string & section ) : RegressionAlgorithmKernel ( conf, kernelFunction )
{
	this->verbose = conf->gB(section, "verbose", false );
	this->optimizeParameters = conf->gB(section, "optimize_parameters", true );
	this->maxIterations = conf->gI( section, "optimization_maxiterations", 30 );
	string optimizationMethod_s = conf->gS(section, "optimization_method", "rasmussen" );

	if ( optimizationMethod_s == "rasmussen" ) 
		optimizationMethod = OPTIMIZATION_METHOD_RASMUSSEN;
	else if ( optimizationMethod_s == "trustregion" )
		optimizationMethod = OPTIMIZATION_METHOD_TRUSTREGION;
	else
		fthrow(Exception, "Optimization method " << optimizationMethod_s << " is unknown.");

	ParameterizedKernel *pkernelFunction = dynamic_cast< ParameterizedKernel * > ( kernelFunction );
	if ( optimizeParameters && (pkernelFunction == NULL) )
	{
		cerr << "RegGaussianProcess: Unable to optimize hyperparameters with no specified kernel function" << endl;
		cerr << "RegGaussianProcess: Switching to non-optimization mode" << endl;
		optimizeParameters = false;
	}

	bool approximateTraceTerm = conf->gB(section, "approximate_trace_term", false);
	if ( approximateTraceTerm ) 
		traceApproximation = new TraceApproximation ( conf, section );
	else
		traceApproximation = NULL;

	useLooParameters = conf->gB(section, "use_loo_parameters", true );

	modelselcrit = NULL;
	if ( useLooParameters ) {
		string modelselcrit_s = conf->gS("section", "loo_crit", "loo_pred_prob" );
		modelselcrit = GenericGPModelSelection::selectModelSel ( conf, modelselcrit_s );
	}

}

RegGaussianProcess::RegGaussianProcess ( const RegGaussianProcess & src ) : 
	RegressionAlgorithmKernel ( src )
{
	kInvY = src.kInvY;
	verbose = src.verbose;
	optimizeParameters = src.optimizeParameters;
	optimizationMethod = src.optimizationMethod;
	traceApproximation = src.traceApproximation;
	modelselcrit = src.modelselcrit;
}

RegGaussianProcess::~RegGaussianProcess()
{
	if ( traceApproximation != NULL )
		delete traceApproximation;
	if ( modelselcrit != NULL )
		delete modelselcrit;

}

void RegGaussianProcess::teach ( KernelData *kernelData, const NICE::Vector & y )
{
	if ( optimizeParameters ) 
	{
		if ( (kernelFunction != NULL) )
		{
			ParameterizedKernel *kernelPara = dynamic_cast< ParameterizedKernel * > ( kernelFunction );
			if ( kernelPara == NULL ) {
				fthrow(Exception, "RegGaussianProcess: you have to specify a parameterized kernel !");
			}
			GPRegressionOptimizationProblem gpopt ( kernelData, y, kernelPara, 
				verbose, modelselcrit, traceApproximation );
			cout << "RegGaussianProcess: Hyperparameter optimization ..." << endl;

			if ( optimizationMethod == OPTIMIZATION_METHOD_TRUSTREGION )
			{
				if ( verbose ) 
					cerr << "RegGaussianProcess: using trust region optimizer" << endl;
				FirstOrderTrustRegion *optimizer = new FirstOrderTrustRegion(0.1 /*typical gradient*/, verbose);
				optimizer->setEpsilonG ( 0.01 );
				optimizer->setMaxIterations ( maxIterations );
				optimizer->optimizeFirst ( gpopt );
				delete optimizer;

			} else if ( optimizationMethod == OPTIMIZATION_METHOD_RASMUSSEN ) {
				if ( verbose ) 
					cerr << "RegGaussianProcess: using conjugate gradient optimizer" << endl;

				FirstOrderRasmussen *optimizer = new FirstOrderRasmussen();
				optimizer->setEpsilonG ( 0.01 );
				optimizer->setMaxIterations ( -maxIterations );
				optimizer->optimizeFirst ( gpopt );
				delete optimizer;
			} else {
				fthrow(Exception, "Unknown optimization method " << optimizationMethod );
			}
			
			cout << "RegGaussianProcess: Hyperparameter optimization ...done" << endl;
		
			if ( useLooParameters ) 
			{
				cerr << "RegGaussianProcess: using best loo parameters" << endl;
				gpopt.useLooParameters();
			}

			gpopt.update();

			Vector parameters;
			kernelPara->getParameters ( parameters );
			cout << "RegGaussianProcess: Optimization finished: " << parameters << endl << endl;
		} else {
			fthrow(Exception, "KCGPRegression: you have to specify a kernel function !" );
		}
	} else {

		if ( !kernelData->hasCholeskyFactorization() )
			kernelData->updateCholeskyFactorization();
	}

	this->y.resize ( y.size() );
	this->y = y;
	this->kInvY.resize ( this->y.size() );

	kernelData->computeInverseKernelMultiply ( y, this->kInvY );
}

double RegGaussianProcess::predictKernel ( const NICE::Vector & kernelVector, double kernelSelf )
{
	// kernelSelf is not needed for the regression type of GP

	if ( kernelVector.size() != this->y.size() ) 
		fthrow(Exception, "RegGaussianProcess::predictKernel: size of kernel value vector " << 
			kernelVector.size() << " does not match number of training points " << this->y.size() );

	double yEstimate = kernelVector.scalarProduct ( kInvY );

	return yEstimate;
}

	
RegGaussianProcess *RegGaussianProcess::clone(void) const
{
	return new RegGaussianProcess ( *this );
}

void RegGaussianProcess::store(std::ostream& ofs, int type) const
{
	ofs << kInvY << endl;
	ofs << this->y << endl;
}

void RegGaussianProcess::restore(std::istream& ifs, int type)
{
	ifs >> kInvY;
	ifs >> this->y;
}