/** 
* @file KCNullSpace.cpp
* @brief Classification and novelty detection with kernel null space methods (Kernel Null Foley Sammon Transform - KNFST)
* @author Paul Bodesheim
* @date 26/11/2012

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

#include "core/vector/Algorithms.h"

#include "KCNullSpace.h"
#include <limits>

#undef DEBUG

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

KCNullSpace::KCNullSpace( const Config *conf, Kernel *kernelFunction, const string & section )
	: KernelClassifier ( conf, kernelFunction )
{
	this->maxClassNo = 0;
	this->verbose = conf->gB( section, "verbose", false );

}

KCNullSpace::KCNullSpace( const KCNullSpace &vcova ): KernelClassifier(vcova)
{

	verbose = vcova.verbose;
        dimNullSpace = vcova.dimNullSpace;
        oneClassSetting = vcova.oneClassSetting;
        
        trainingSetStatistic.clear();
        for ( std::map<int,int>::const_iterator it = vcova.trainingSetStatistic.begin(); it != vcova.trainingSetStatistic.end(); it++ )
        {
          trainingSetStatistic.insert(pair<int,int>( (*it).first,(*it).second ));
        }
	
        nullProjectionDirections.resize( vcova.nullProjectionDirections.rows(),vcova.nullProjectionDirections.cols() );
        nullProjectionDirections.set( 0.0 );
	for(int i = 0; i < (int)vcova.nullProjectionDirections.rows(); i++)
	{
          for(int j = 0; j < (int)vcova.nullProjectionDirections.cols(); j++)
          {
            nullProjectionDirections(i,j) = vcova.nullProjectionDirections(i,j);
          }          
	}

        targetPoints.clear();
        for( std::map<int,NICE::Vector>::const_iterator it = vcova.targetPoints.begin(); it != vcova.targetPoints.end(); it++ )
        {
          targetPoints.insert(pair<int,NICE::Vector>( (*it).first,(*it).second ));
        }
        
        eigenBasis.resize( vcova.eigenBasis.rows(),vcova.eigenBasis.cols() );
        eigenBasis.set( 0.0 );
        for(int i = 0; i < (int)vcova.eigenBasis.rows(); i++)
        {
          for(int j = 0; j < (int)vcova.eigenBasis.cols(); j++)
          {
                eigenBasis(i,j) = vcova.eigenBasis(i,j);
          }
        }
}

KCNullSpace::~KCNullSpace()
{

}

void KCNullSpace::teach ( KernelData *kernelData, const NICE::Vector & y )
{
        NICE::Vector labels(y);
        int maxLabel( (int)labels.Max());
        
        /** check if we are in a one-class setting */
        int minLabel( (int)labels.Min());
        if (maxLabel == minLabel)
        {
          oneClassSetting = true;
          maxClassNo = 1;
          if (verbose)
            std::cerr << "KCNullSpace::teach: one-class setting" << std::endl;
          /** one-class setting: add a row and a column of zeros to the kernel matrix representing dot products with origin in kernel feature space*/
          kernelData->increase_size_by_One();
          kernelData->getKernelMatrix() (kernelData->getKernelMatrix().rows()-1, kernelData->getKernelMatrix().cols()-1) = 0.0;
          labels.append(minLabel+1);
          computeTrainingSetStatistic(labels);
          
        }
        else
        {
          oneClassSetting = false;
          if (verbose)
            std::cerr << "KCNullSpace::teach: multi-class setting" << std::endl;          
          computeTrainingSetStatistic(labels);
          maxClassNo = trainingSetStatistic.size()-1; // number of classes, start counting at 0
        }
         
        computeNullProjectionDirections(kernelData,labels);
 
        computeTargetPoints(kernelData,labels);
        
        if (oneClassSetting)
        {
          /** remove the value that corresponds to the dot product with the origin in the kernel feature space, since this is always equal to zero */
          nullProjectionDirections.deleteRow(nullProjectionDirections.rows()-1);
        }
}

std::map<int,int> * KCNullSpace::getTrainingSetStatistic()
{
  return &trainingSetStatistic; 
}

NICE::Matrix KCNullSpace::getNullProjectionDirections()
{
  return nullProjectionDirections;
}

std::map<int,NICE::Vector> *  KCNullSpace::getTargetPoints()
{
  return &targetPoints;
}

int  KCNullSpace::getNullSpaceDimension()
{
  return dimNullSpace;
}

bool  KCNullSpace::isOneClass()
{
  return oneClassSetting;
}

void KCNullSpace::computeTrainingSetStatistic(const NICE::Vector & y)
{
        trainingSetStatistic.clear();
        std::map<int,int>::iterator it;
        for ( uint i = 0 ; i < y.size(); i++ )
        {
          it = trainingSetStatistic.find ( (int)y[i] );
          if ( it == trainingSetStatistic.end() )
          {
            trainingSetStatistic.insert(it, pair<int,int>( (int)y[i],1 ));
          }
          else
          {
            it->second += 1;
          }
        }  
}

void KCNullSpace::computeBasisUsingKernelPCA(const KernelData *kernelData)
{
        if (verbose)
          std::cerr << "KCNullSpace::computeBasisUsingKernelPCA: compute kernel PCA basis..." << std::endl;
        
        NICE::Matrix K (kernelData->getKernelMatrix());
        
        /** let K represent dot products of zero mean data in kernel feature space */
        centerKernelMatrix(K);

        /** get eigenvectors and eigenvalues (descreasing order) of centered kernel matrix*/
        NICE::Matrix eigenVectors(K.rows(), K.cols(), 0.0);
        NICE::Vector eigenValues(K.rows(), 0.0);
        K.addIdentity(1.0);
        eigenvectorvalues(K, eigenVectors, eigenValues);
        eigenValues -= 1.0;
//         NICE::GMCovariance gm(&K);
//         NICE::EigValuesTRLAN ev;
//         ev.getEigenvalues(gm, eigenValues, eigenVectors, K.rows());
        
        /** only use eigenvectors of non-zero eigenvalues */
        int j(0);
        for (size_t i=0; i<K.rows(); i++)
        {
          if ( eigenValues(i) < 1e-12 )
          {
            eigenVectors.deleteCol(i);
          } 
          else
          {
            j++;
          }
        }
        eigenValues.resize(j);
        
        /** scale eigenvectors with eigenvalues */
        double scale(0.0);
        for (size_t c=0; c<eigenVectors.cols(); c++)
        {
          scale = 1.0/sqrt(eigenValues(c));
          for (size_t r=0; r<eigenVectors.rows(); r++)
            eigenVectors(r,c) *= scale;
        }
        
        eigenBasis = eigenVectors;
        
        if (verbose)
          std::cerr << "KCNullSpace::computeBasisUsingKernelPCA: computation done" << std::endl;        
    
}
    
void KCNullSpace::centerKernelMatrix(NICE::Matrix & kernelMatrix)
{
        NICE::Matrix onesK (kernelMatrix.rows(), kernelMatrix.cols(), 1.0/kernelMatrix.rows());
        kernelMatrix = kernelMatrix - onesK*kernelMatrix - kernelMatrix*onesK + onesK*kernelMatrix*onesK;
}
    
void KCNullSpace::computeNullProjectionDirections ( const KernelData *kernelData, const NICE::Vector & y )
{
         
        if (verbose)
          std::cerr << "KCNullSpace::computeNullProjectionDirections: compute null projection directions..." << std::endl; 
        
        /** obtain Kernel PCA basis */
        computeBasisUsingKernelPCA(kernelData);      
        
        /** set matrix IM=(I-M) with I being the unit matrix and M being a matrix with all entries equal to 1/n where n is the number of training samples */
        NICE::Matrix IM (y.size(),y.size(),-1.0/y.size());
        IM.addIdentity(1.0);
        
        /** compute matrix IL=(I-L) with I being the unit matrix and L being a block matrix where in each block of class samples each entry is equal to 1/numClassSamples */
        NICE::Matrix IL (y.size(),y.size(),0.0);
        IL.addIdentity(1.0);        
        for (size_t c=0; c<IL.cols(); c++)
        {
          /** if sample with index r is in the same class as sample with index c, then insert the value 1/numClassSamples */
          for (size_t r=0; r<IL.rows(); r++)
          { 
            if ( y(r) == y(c) )
            {
              IL(r,c) -= 1.0/trainingSetStatistic[(int)y(r)];;
            }
          }
        }
        
        /** compute Matrix H = ((I-M)*basisvecs)^T * kernelMatrix * (I-L) with I being the unit matrix */
        NICE::Matrix H (y.size(),y.size(),0.0);
        IM = IM*eigenBasis;
        H = IM.transpose() * kernelData->getKernelMatrix() * IL;
        
        /** obtain matrix T = H * H^T */
        NICE::Matrix T = H*H.transpose();
        
        /** get eigenvectors and eigenvalues (descreasing order) of T */
        NICE::Matrix eigenVectors(T.rows(), T.cols(), 0.0);
        NICE::Vector eigenValues(T.rows(), 0.0);
        T.addIdentity(1.0);
        eigenvectorvalues(T, eigenVectors, eigenValues);
        eigenValues -= 1.0;
//         NICE::GMCovariance gm(&T);
//         NICE::EigValuesTRLAN ev;
//         ev.getEigenvalues(gm, eigenValues, eigenVectors, T.rows());
        
        if (verbose)
          std::cerr << "T: " << T.rows() << " x " << T.rows() << " nan: " << T.containsNaN()<< std::endl;
        
        /** only use eigenvectors of zero eigenvalues (null space!!!) but at least one eigenvector according to the smallest eigenvalue (therefore start at index i=T.rows()-2)*/
        for (int i=T.rows()-2; i>=0; i--)
        {
          if ( eigenValues(i) > 1e-12 )
          {
            eigenVectors.deleteCol(i);
          } 
        }
        
        /** compute null projection directions */
        nullProjectionDirections = IM*eigenVectors;
        dimNullSpace = nullProjectionDirections.cols();
        
        if (verbose)
          std::cerr << "KCNullSpace::computeNullProjectionDirections: computation done" << std::endl;
}

void KCNullSpace::computeTargetPoints ( const KernelData *kernelData, const NICE::Vector & y )
{
  if (verbose)
    std::cerr << "KCNullSpace::computeTargetPoints: compute target points..." << std::endl; 
          
  targetPoints.clear();
  NICE::Vector targetPoint (dimNullSpace, 0.0);
  int classLabel(0);
  
  if (oneClassSetting)
  {
   
    std::map<int,int>::iterator it;
    it = trainingSetStatistic.begin();
    classLabel = (*it).first;
    for (size_t i=0; i<y.size(); i++)
    {
      if ( classLabel == y(i) )
      {
        /** project each training sample in the null space */
        targetPoint += nullProjectionDirections.transpose() * kernelData->getKernelMatrix().getColumn(i);
      }
    }
    /** average the projections of all class samples due to numerical noise */
    targetPoint /= (*it).second;
    
    /** we only have one target point in an one-class setting */
    targetPoints.insert(pair<int,NICE::Vector>( classLabel,targetPoint));
    
  }
  else
  {
    
    /** create averaging vectors for each class, necessary to compute target points at the end of this method */
    std::map<int,NICE::Vector> averagingVectors;
    averagingVectors.clear();
    NICE::Vector averagingVector(y.size(),0.0);
    
    /** insert one averaging vector  for each class */
    int numClassSamples(0);
    for ( std::map<int,int>::iterator it = trainingSetStatistic.begin(); it != trainingSetStatistic.end(); it++ )
    {
      
      /** create current averaging vector */
      classLabel = (*it).first;
      numClassSamples = (*it).second;
      for ( size_t j = 0; j < y.size(); j++ )
      {
        if ( classLabel == y[j]  )
        {
          averagingVector(j) = 1.0/numClassSamples;
        }
        else
        {
          averagingVector(j) = 0.0;
        }
      }
      
      /** insert averaging vector for current class*/
      averagingVectors.insert(pair<int,NICE::Vector>( classLabel,averagingVector));
    }
    
    /** compute target points using previously created averaging vectors: average for each class the projections of the class samples in the null space */
    for ( std::map<int,NICE::Vector>::iterator it = averagingVectors.begin(); it != averagingVectors.end(); it++ )
    {
       targetPoint = nullProjectionDirections.transpose() * kernelData->getKernelMatrix() * (*it).second;
       targetPoints.insert(pair<int,NICE::Vector>( (*it).first,targetPoint));
    }    
    
  }
  
  if (verbose)
    std::cerr << "KCNullSpace::computeTargetPoints: computation done" << std::endl;
  
}

ClassificationResult KCNullSpace::classifyKernel ( const NICE::Vector & kernelVector, double kernelSelf ) const
{
	if ( targetPoints.size() <= 0 ) 
		fthrow(Exception, "The classifier was not trained with training data (use teach(...))");
        
        if (oneClassSetting)
        {
          return noveltyDetection ( kernelVector, kernelSelf );
        }  
        
        NICE::Vector projection(dimNullSpace,0.0);
        projection = nullProjectionDirections.transpose() * kernelVector;
        
	FullVector scores ( maxClassNo+1 );
	scores.set(- std::numeric_limits<double>::max());
        
        int iter(0);     
        for ( std::map<int,NICE::Vector>::const_iterator it = targetPoints.begin(); it != targetPoints.end(); it++ )
        {
          scores[iter] = 1.0-(it->second - projection).normL2();
          iter++;         
        }
        
	ClassificationResult r ( scores.maxElement(), scores );
	return r;
}

ClassificationResult KCNullSpace::noveltyDetection ( const NICE::Vector & kernelVector, double kernelSelf ) const
{
        if ( targetPoints.size() <= 0 ) 
                fthrow(Exception, "The classifier was not trained with training data (use teach(...))");
               
        NICE::Vector projection(dimNullSpace,0.0);
        projection = nullProjectionDirections.transpose() * kernelVector;
        
        FullVector scores ( 2 );
        scores.set(- std::numeric_limits<double>::max());
        double tmp_score(0.0);
        
        for ( std::map<int,NICE::Vector>::const_iterator it = targetPoints.begin(); it != targetPoints.end(); it++ )
        {
          tmp_score = 1.0-(it->second - projection).normL2();
          if (tmp_score > scores[1])
          {
            scores[1] = tmp_score;
          }
        }
        
        ClassificationResult r ( scores.maxElement(), scores );
        return r;
}

KCNullSpace* KCNullSpace::clone(void) const
{
	KCNullSpace *classifier = new KCNullSpace( *this );
	return classifier;
}

void KCNullSpace::clear()
{
	//nothing to clear
}

void KCNullSpace::restore(std::istream& ifs, int type)
{
	ifs.precision (numeric_limits<double>::digits10 + 1);
	ifs >> maxClassNo;
	ifs >> oneClassSetting;
	ifs >> dimNullSpace;

        ifs >> eigenBasis;
        
        int k(0);
        int classLabel(0);
        int numClassSamples(0);
        ifs >> k;
        trainingSetStatistic.clear();
        for (int i=0; i<k; i++)
        {
          ifs >> classLabel;
          ifs >> numClassSamples;
          trainingSetStatistic.insert( pair<int,int>(classLabel,numClassSamples) );
        }
               
        ifs >> nullProjectionDirections;
        
        NICE::Vector targetPoint;
        ifs >> k;
        for (int i=0; i<k; i++)
        {
          ifs >> classLabel;
          ifs >> targetPoint;
          targetPoints.insert( pair<int,NICE::Vector>(classLabel,targetPoint) );
        }
        
	KernelClassifier::restore(ifs,type);
}

void KCNullSpace::store(std::ostream& ofs, int type) const
{
	ofs.precision (numeric_limits<double>::digits10 + 1);
	
	ofs << maxClassNo << endl;
	ofs << oneClassSetting << endl;
	ofs << dimNullSpace << endl;
        
        ofs << eigenBasis << endl;
        
        ofs << trainingSetStatistic.size() << endl;
        for (std::map<int,int>::const_iterator it = trainingSetStatistic.begin() ; it != trainingSetStatistic.end(); it++)
        {
          ofs << (*it).first << endl;
          ofs << (*it).second << endl;
        }
        
        ofs << nullProjectionDirections << endl;
        
        ofs << targetPoints.size() << endl;
        for (std::map<int,NICE::Vector>::const_iterator it = targetPoints.begin() ; it != targetPoints.end(); it++)
        {
          
          ofs << (*it).first << endl;
          ofs << (*it).second << endl;
        }
	
	KernelClassifier::store(ofs,type);
}