/** 
* @file KernelLinCombNormalized.cpp
* @brief Interface for the popular exponential mercer kernels
* @author Erik Rodner
* @date 10/24/2007

*/

#include <iostream>

#include <math.h>
#include "KernelLinCombNormalized.h"

using namespace OBJREC;

using namespace std;
using namespace NICE;



KernelLinCombNormalized::KernelLinCombNormalized( uint num, bool _exponential )
	: alpha ( num ), exponential(_exponential)
{
	for ( uint i = 0 ; i < num ; i++ )
		alpha[i] = - log((double)num);
}

KernelLinCombNormalized::KernelLinCombNormalized ( const KernelLinCombNormalized & src )
{
	alpha.resize ( src.alpha.size() );
	alpha = src.alpha;
	exponential = src.exponential;
}

KernelLinCombNormalized::~KernelLinCombNormalized()
{
}

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

	
double KernelLinCombNormalized::K (const NICE::Vector & x, const NICE::Vector & y) const
{
	fthrow(Exception, "Unable to calculate kernel values directly.");
}

void KernelLinCombNormalized::updateKernelData ( KernelData *kernelData ) const
{
	NICE::Matrix & kernelMatrix = kernelData->getKernelMatrix();

	double norm = 0.0;
	for ( uint k = 0 ; k < alpha.size(); k++ )
		norm += exp(alpha[k]);

	for ( uint k = 0 ; k < alpha.size(); k++ )
	{
		const NICE::Matrix & cachedMatrix = kernelData->getCachedMatrix ( k );
			
		if ( k == 0 ) {
			kernelMatrix.resize ( cachedMatrix.rows(), cachedMatrix.cols() );
			kernelMatrix.set(0.0);
		}
		kernelMatrix += exp(alpha[k])/norm * cachedMatrix;
	}

	if ( exponential ) {
		for ( uint i = 0 ; i < kernelMatrix.rows(); i++ )
			for ( uint j = 0 ; j < kernelMatrix.cols(); j++ )
				kernelMatrix(i,j) = exp(- kernelMatrix(i,j));
	}
}

void KernelLinCombNormalized::getKernelJacobi ( size_t parameter, const NICE::Vector & parameters, const KernelData *kernelData, NICE::Matrix & jacobiMatrix ) const
{
	// double _gammasq = exp(2*parameters[0]);
	const NICE::Matrix & kernelMatrix = kernelData->getKernelMatrix();

	if ( parameter >= alpha.size() )
		fthrow(Exception, "Parameter index exceeds parameter vector size" );

	jacobiMatrix.resize ( kernelMatrix.rows(), kernelMatrix.cols() );

#if 0
	if ( exponential ) {
		for ( uint i = 0 ; i < jacobiMatrix.rows(); i++ )
			for ( uint j = 0 ; j < jacobiMatrix.cols(); j++ )
				jacobiMatrix(i,j) = - kernelMatrixParameter(i,j) * exp(alpha[parameter]) * kernelMatrix(i,j);
	} else {
		jacobiMatrix = exp(alpha[parameter]) * kernelMatrixParameter;
	}
#endif
	// derivation of \sum_k exp(alpha_k) / (sum_j exp(alpha_j)) * K_k = \sum_k exp(\alpha_k) / norm * K_k
	//
	// \sum_k (k==parameter) ...
	
	double norm = 0.0;
	for ( uint k = 0 ; k < alpha.size(); k++ )
		norm += exp(alpha[k]);

	Matrix tmp ( kernelMatrix.rows(), kernelMatrix.cols() );
	tmp.set(0.0);
	for ( uint k = 0 ; k < alpha.size(); k++ )
	{
		const Matrix & kernelMatrixSingle = kernelData->getCachedMatrix( k );
		for ( uint i = 0 ; i < tmp.rows(); i++ )
			for ( uint j = 0 ; j < tmp.cols(); j++ )
				if ( k == parameter ) {
					tmp(i,j) += (exp(alpha[k])/norm) * ( 1 - exp(alpha[k])/norm ) * kernelMatrixSingle(i,j);
				} else {
					tmp(i,j) -= (exp(alpha[k])/norm) * (exp(alpha[parameter])/norm) * kernelMatrixSingle(i,j);
				}
	}

	if ( exponential ) {
		for ( uint i = 0 ; i < jacobiMatrix.rows(); i++ )
			for ( uint j = 0 ; j < jacobiMatrix.cols(); j++ )
				jacobiMatrix(i,j) = - tmp(i,j) * kernelMatrix(i,j);
	} else {
		jacobiMatrix = tmp;
	}


}

void KernelLinCombNormalized::setParameters( const NICE::Vector & newParameters ) 
{
	alpha = newParameters;
}

void KernelLinCombNormalized::getParameters( NICE::Vector & newParameters ) const
{
	newParameters.resize ( alpha.size() );
	newParameters = alpha;
}