/**
 * @file PCA.cpp
 * @brief Computation of a PCA by Eigenvalue Decomposition, Karhunen Loewe Transform
 * @author Michael Koch
 * @date 5/27/2008

 */
/** \example testPCA.cpp
 * This is an example of how to use the Principal Component Analysis (PCA) class.
 */

// #define DEBUG
#undef DEBUG

#include <iostream>
#include <math.h>
#include <vector>
#include <map>

#include "core/algebra/GenericMatrix.h"
#include "core/algebra/EigValues.h"
#include "core/algebra/EigValuesTRLAN.h"

// #include "vislearning/fourier/FourierLibrary.h"
#include "PCA.h"
#include "vislearning/baselib/Gnuplot.h"

using namespace OBJREC;

using namespace std;
using namespace NICE;

PCA::PCA(uint dim, uint maxiteration, double mindelta)
{
	init(dim, maxiteration, mindelta);
}

PCA::PCA(void)
{
	init();
}

PCA::~PCA()
{
}

void PCA::init(uint dim, uint maxiteration, double mindelta)
{
	this->targetDimension = dim;
	this->maxiteration = maxiteration;
	this->mindelta = mindelta;
}

void PCA::restore(istream & is, int format)
{

	is >> basis;
	is >> normalization;
	is >> mean;
	is >> targetDimension;
}

void PCA::store(ostream & os, int format) const
{
	os << basis << normalization << mean << targetDimension;
}

void PCA::clear()
{

}

void PCA::calculateBasis(const NICE::Matrix &features,
		const uint targetDimension, const uint mode)
{
	calculateBasis(features, targetDimension, false);
}

void PCA::calculateMean ( const NICE::Matrix &features, NICE::Vector & mean )
{
  // data vectors are put row-wise in the matrix
  mean.resize(features.cols());
  for ( uint i = 0 ; i < features.rows(); i++ )
    mean = mean + features.getRow(i);
}

void PCA::calculateBasis(const NICE::Matrix &features,
		const uint targetDimension, const bool adaptive,
		const double targetRatio)
{
	this->targetDimension = targetDimension;
  // dimension of the feature vectors
	uint srcFeatureSize = features.cols();
	uint mindimension = std::min(this->targetDimension, srcFeatureSize);

	NICE::Matrix eigenvectors;

	NICE::Vector eigenvalue;
  
  calculateMean(features, mean);
#ifdef NICE_USELIB_TRLAN
  EigValues *eig = new EigValuesTRLAN();//fast lanczos TRLAN
#else
  EigValues *eig = new EVArnoldi();//Arnoldi for (srcFeatureSize<n)
#endif
  NICE::Matrix features_transpose = features.transpose();
  GMCovariance C(&features_transpose);
  if (adaptive)
  {
    eig->getEigenvalues(C, eigenvalue, eigenvectors, srcFeatureSize);
  }
  else
  {
    eig->getEigenvalues(C, eigenvalue, eigenvectors, mindimension);
  }

#ifdef DEBUG
	fprintf(stderr, "Eigenvalue Decomposition ready \n");
	cerr << eigenvectors << endl;
	cerr << eigenvalue << endl;

	//sort values
	fprintf(stderr, "Eigenvector-Rows:%i Eigenvector-Cols:%i\n", (int)eigenvectors.rows(), (int)eigenvectors.cols());
#endif

	multimap<double, NICE::Vector> map;
	double sumeigen = 0.0;

	NICE::Vector ratio(eigenvectors.cols());
	for (uint i = 0; i < eigenvectors.cols(); i++)//every eigenvector
	{
		NICE::Vector eigenvector(srcFeatureSize);
		for (uint k = 0; k < srcFeatureSize; k++)
		{
			eigenvector[k] = eigenvectors(k, i);
		}
		map.insert(pair<double, NICE::Vector> (eigenvalue[i], eigenvector));
		sumeigen += eigenvalue[i];
	}

	//compute basis size


	if (adaptive)
	{ //compute target dimension
		uint dimensioncount = 0;
		double addedratio = 0.0;
		multimap<double, NICE::Vector>::reverse_iterator it = map.rbegin();
		while (addedratio <= targetRatio && it != map.rend())
		{
			//calc ratio
			ratio[dimensioncount] = (*it).first / sumeigen;
			addedratio += ratio[dimensioncount];
			dimensioncount++;
			it++;
		}
		this->targetDimension = dimensioncount;
	}

	mindimension = std::min(this->targetDimension, srcFeatureSize);
	this->targetDimension = mindimension;
	basis = NICE::Matrix(srcFeatureSize, mindimension);
	//get sorted values
	uint count = 0;
	multimap<double, NICE::Vector>::reverse_iterator it = map.rbegin();
	while (count < this->targetDimension && it != map.rend())
	{
		NICE::Vector eigenvector = (*it).second;
		//put eigenvector into column
		for (uint k = 0; k < srcFeatureSize; k++)
		{
			basis(k, count) = eigenvector[k];
		}
		//calc ratio
		ratio[count] = (*it).first / sumeigen;

		count++;
		it++;
	}
	//normalization matrix / modify variance to 1 for all eigenvectors
	normalization = NICE::Matrix(mindimension, mindimension, 0);
	for (uint k = 0; k < mindimension; k++)
	{
		normalization(k, k) = 1.0 / sqrt(eigenvalue[k]);
	}

#ifdef DEBUG
	cout << "Eigenvalue-absolute:" << eigenvalue << endl;
	cout << "Eigenvalue-ratio:" << ratio << endl;
#endif

}

NICE::Vector PCA::getFeatureVector(const NICE::Vector &data,
		const bool normalize)
{
	//free data of mean
	if (normalize)
	{

		NICE::Vector meanfree(data);
		meanfree -= mean;
		//Y=W^t * B^T
		if(normbasis.rows() == 0)
			normbasis.multiply(normalization, basis, false, true);
		NICE::Vector tmp;
		tmp.multiply(normbasis, meanfree);
		return tmp;
	}
	else
	{
		NICE::Vector meanfree(data);
		meanfree -= mean;
		//Y=W^t * B^T
		NICE::Vector tmp;
		tmp.multiply(basis, meanfree, true);
		return tmp;
	}
}

NICE::Vector PCA::getMean()
{
	return mean;
}

NICE::Matrix PCA::getBasis()
{
	return basis;
}

int PCA::getTargetDim()
{
	return targetDimension;
}