/** 
* @file KCNullSpaceNovelty.cpp
* @brief 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 "KCNullSpaceNovelty.h"
#include <limits>

#undef DEBUG

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

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

}

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

	verbose = vcova.verbose;
        dimNullSpace = vcova.dimNullSpace;
        oneClassSetting = vcova.oneClassSetting;
        
        trainingSetStatistic.clear();
        std::map<int,int>::iterator it;
        for ( it = ( (std::map<int,int>)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(int i = 0; i < (int)vcova.targetPoints.size(); i++)
        {
                targetPoints.push_back( NICE::Vector(vcova.targetPoints[i]) );
        }
        
        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);
          }
        }
}

KCNullSpaceNovelty::~KCNullSpaceNovelty()
{

}

void KCNullSpaceNovelty::teach ( KernelData *kernelData, const NICE::Vector & y )
{
        NICE::Vector labels(y);
	maxClassNo = (int)labels.Max();
        
        /** check if we are in a one-class setting */
        int minClassNo = (int)labels.Min();
        if (maxClassNo == minClassNo)
        {
          oneClassSetting = true;
          if (verbose)
            std::cerr << "KCNullSpaceNovelty::teach: one-class setting" << std::endl;
          computeTrainingSetStatistic(labels);
          /** 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(minClassNo+1);
        }
        else
        {
          oneClassSetting = false;
          if (verbose)
            std::cerr << "KCNullSpaceNovelty::teach: multi-class setting" << std::endl;          
          computeTrainingSetStatistic(labels);
        }

        if (verbose)
          std::cerr << "KCNullSpaceNovelty::teach: compute null projection directions..." << std::endl;   
        computeNullProjectionDirections(kernelData,labels);
        if (verbose)
          std::cerr << "KCNullSpaceNovelty::teach: compute target points..." << std::endl;  
        computeTargetPoints(kernelData,labels);
}

std::map<int,int> KCNullSpaceNovelty::getTrainingSetStatistic()
{
  return trainingSetStatistic; 
}

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

NICE::VVector  KCNullSpaceNovelty::getTargetPoints()
{
  return targetPoints;
}

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

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

void KCNullSpaceNovelty::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 KCNullSpaceNovelty::computeBasisUsingKernelPCA(const KernelData *kernelData)
{
        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);
        eigenvectorvalues(K, eigenVectors, eigenValues);

        /** 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;
    
}
    
void KCNullSpaceNovelty::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 KCNullSpaceNovelty::computeNullProjectionDirections ( const KernelData *kernelData, const NICE::Vector & y )
{
         
        /** obtain Kernel PCA basis */
        if (verbose)
          std::cerr << "KCNullSpaceNovelty::computeNullProjectionDirections: compute kernel PCA basis..." << std::endl;
        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; 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);
        eigenvectorvalues(T, eigenVectors, eigenValues);
        
        /** 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 << "KCNullSpaceNovelty::computeNullProjectionDirections: computation done" << std::endl;
}

void KCNullSpaceNovelty::computeTargetPoints ( const KernelData *kernelData, const NICE::Vector & y )
{
  
  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.push_back(targetPoint);
    
  }
  else
  {
    
    /** create averaging vectors for each class, necessary to compute target points at the end of this method */
    std::vector<NICE::Vector> averagingVectors;
    averagingVectors.clear();
    NICE::Vector averagingVector(y.size(),0.0);
    
    /** insert one averaging vector and one target point vector for each class */
    std::map<int,int>::iterator it;
    int numClassSamples(0);
    for ( 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.push_back(averagingVector);
      /** insert a null vector for the target point of the current class */
      targetPoints.push_back(targetPoint);
    }
    
    /** compute target points using previously created averaging vectors: average for each class the projections of the class samples in the null space */
    for ( size_t i = 0 ; i < targetPoints.size(); i++ )
    {
      targetPoints[i] = nullProjectionDirections.transpose() * kernelData->getKernelMatrix() * averagingVectors[i];
    }    
    
  }
  
  if (verbose)
    std::cerr << "KCNullSpaceNovelty::computeTargetPoints: computation done" << std::endl;
  
}

ClassificationResult KCNullSpaceNovelty::classifyKernel ( 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 ( trainingSetStatistic.size() );
	scores.set(0);
        
        std::map<int,int>::iterator it;
        int iter(0);
        for ( it = ( (std::map<int,int>)trainingSetStatistic ).begin(); it != trainingSetStatistic.end(); it++ )
        {
          scores[iter] = -(targetPoints[iter] - projection).normL2();
          iter++;         
        }
        
	ClassificationResult r ( scores.maxElement(), scores );

	return r;
}

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

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

void KCNullSpaceNovelty::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 >> targetPoint;
          targetPoints.push_back(targetPoint);
        }
        
	KernelClassifier::restore(ifs,type);
}

void KCNullSpaceNovelty::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;
        std::map<int,int>::iterator it;
        for (it = ( (std::map<int,int>)trainingSetStatistic ).begin() ; it != trainingSetStatistic.end(); it++)
        {
          ofs << (*it).first << endl;
          ofs << (*it).second << endl;
        }
        
        ofs << nullProjectionDirections << endl;
        
        ofs << targetPoints.size() << endl;
        for (size_t k=0; k<targetPoints.size(); k++)
        {
          ofs << targetPoints[k] << endl;
        }
	
	KernelClassifier::store(ofs,type);
}