/**
 * @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() );
  mean.set(0);
  for ( uint i = 0 ; i < features.rows(); i++ )
    mean = mean + features.getRow ( i );
  mean /= features.rows();
}

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;
  if(mindimension < 12)
    eig = new EVArnoldi();//Arnoldi for (srcFeatureSize<n)
  else
    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;
}