/** 
* @file SpectralCluster.cpp
* @brief spectral clustering by kmeans-clustering of eigenvectors
* @author Erik Rodner
* @date 11/13/2007

*/

#include <iostream>

#include <core/vector/Distance.h>
#include <core/vector/Eigen.h>
#include <map>
#include "vislearning/math/cluster/SpectralCluster.h"

using namespace OBJREC;

using namespace std;
using namespace NICE;

SpectralCluster::SpectralCluster ( int _noClusters, double alpha ) : noClusters(_noClusters), kmeans(_noClusters)
{
    this->alpha = alpha;
}

SpectralCluster::SpectralCluster( const NICE::Config *conf, const std::string & _section) : kmeans(conf)
{  
  this->noClusters = conf->gI( _section, "noClusters", 20);
  this->alpha = conf->gD( _section, "alpha", 1.0);
}

SpectralCluster::~SpectralCluster()
{
}

void SpectralCluster::getSimilarityMatrix ( const VVector & features, 
                                            NICE::Matrix & laplacian,
                                            double alpha )
{
    NICE::Matrix distances ( laplacian );
    double mindist = numeric_limits<double>::max();
    double maxdist = - numeric_limits<double>::max();
    long int count = 0;
    double mean = 0.0;
    double stddev = 0.0;
    NICE::EuclidianDistance<double> distance;

    for ( int i = 0 ; i < (int)features.size() ; i++ )
    {
      const NICE::Vector & xi = features[i];
      for ( int j = i ; j < (int)features.size() ; j++ )
      {
          const NICE::Vector & xj = features[j];
          // double sim = xi * xj;
          double dist = distance.calculate ( xi, xj );
          distances(i, j) = dist;

          if ( dist < mindist )
        mindist = dist;
          if ( dist > maxdist )
        maxdist = dist;

          count++;
          mean += dist;
      }
    }
    mean /= count;

    for ( int i = 0 ; i < (int)features.size() ; i++ )
    {
      for ( int j = i ; j < (int)features.size() ; j++ )
      {
          double d = ( mean - distances(i, j) );
          stddev += d*d;
      }
    }

    stddev /= count;

    double norm = alpha / (2.0 * stddev);

    for ( int i = 0 ; i < (int)features.size() ; i++ )
    {
      for ( int j = i ; j < (int)features.size() ; j++ )
      {
          double sim = exp(-distances(i, j) * norm );
          
          laplacian(i, j) = - sim;
          laplacian(j, i) = - sim;
      }
    }
}

void SpectralCluster::computeLaplacian ( const NICE::VVector & features,
                                         NICE::Matrix & laplacian,
                                         int method, double alpha )
{
    laplacian.set(0.0);
    getSimilarityMatrix(features, laplacian, alpha);
    
    NICE::Vector d ( laplacian.rows() );
    d.set(0.0);
    for ( int i = 0 ; i < (int)laplacian.rows(); i++ )
    {
      for ( int j = 0 ; j < (int)laplacian.cols(); j++ )
      {
          d[i] -=laplacian(i, j);
      }
    }

    // Now we got the negative weight matrix laplacian : -W
    // and the degree matrix : D

        // calculate normalization
    // D^-1 * L
    if ( method == L_RW_NORMALIZED ) 
    {
    // L = D^-1 L_unnormalized = I - D^-1*W
      for ( int i = 0 ; i < (int)laplacian.rows() ; i++ )
      {
        for ( int j = 0 ; j < (int)laplacian.cols() ; j++ )
        {
          laplacian(i, j) *= (1.0/d[i]);
        }
      }

      for ( int i = 0 ; i < (int)laplacian.rows() ; i++ )
      {
          laplacian(i, i) += 1.0;
      }

    }
    else if ( method == L_UNNORMALIZED )
    {
    // unnormalized version
    // L = D - W
      for ( int i = 0 ; i < (int)laplacian.rows() ; i++ )
      {
          laplacian(i, i) += d[i];
      }
    }
}

void SpectralCluster::cluster ( const NICE::VVector & features,
                                NICE::VVector & prototypes,
                                std::vector<double> & weights,
                                std::vector<int>    & assignment )
{
    if ( features.size() <= 0 ) {
      fprintf (stderr, "FATAL ERROR: not enough features vectors provided\n");
      exit(-1);
    }

    const NICE::Vector & x = features[0];
    int dimension = x.size();

    NICE::Matrix laplacian ( features.size(), features.size() );

    computeLaplacian ( features, laplacian, L_RW_NORMALIZED, alpha );
    //computeLaplacian ( features, laplacian, L_UNNORMALIZED, alpha );

    NICE::Matrix eigvect;
    NICE::Vector eigvals;
    try {
      NICE::eigenvectorvalues ( laplacian, eigvect, eigvals );
    }
    catch (...)
    {
      std::cerr << "You should adjust the alpha parameter, or the features are somehow weird" << std::endl;
      std::cerr << "Laplacian = " << laplacian(0, 0, min((int)laplacian.rows()-1, 5), min((int)laplacian.cols()-1, 5)) << std::endl;
      throw Exception ("Laplacian matrix is singular or not well-conditioned!");
    }

    std::map<double, int> eigvals_sorted;
    for ( int i = 0 ; i < (int)eigvals.size(); i++ ) 
    {
      eigvals_sorted.insert ( make_pair( eigvals[i], i ) );
    }

    NICE::VVector spectral_features;

    for ( int i = 0 ; i < (int)eigvect.rows() ; i++ )
    {
      NICE::Vector eigvec_k ( noClusters );
      map<double, int>::const_iterator k = eigvals_sorted.begin();
      for ( int j = 0 ; j < noClusters ; j++ )
      {
        int eigval_index = k->second;
        eigvec_k[j] = eigvect(i, eigval_index ) ;
        k++;
      }

      spectral_features.push_back ( eigvec_k );
    }

    this->kmeans.cluster ( spectral_features, prototypes, weights, assignment ); 

    // recompute prototypes

    for ( int k = 0 ; k < noClusters ; k++ )
    {
      prototypes[k].resize( dimension );
      prototypes[k].set(0);
      weights[k] = 0;
    }

    int j = 0;
    for ( VVector::const_iterator i  = features.begin();
                                  i != features.end();
                                  i++, j++ )
    {
      int k = assignment[j];
      
      NICE::Vector & p = prototypes[k];
      const NICE::Vector & x = *i;
      p += x;
      weights[k]++;
    }

    for ( int k = 0 ; k < noClusters ; k++ )
    {
      NICE::Vector & p = prototypes[k];
      if ( weights[k] <= 0 )
      {
          fprintf (stderr, "FATAL ERROR: spectral clustering produced empty cluster !\n");
          exit(-1);
      }
      p *= ( 1.0 / weights[k] );
      weights[k] = weights[k] / features.size();
    }
}