/** 
* @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;
// refactor-nice.pl: check this substitution
// old: using namespace ice;
using namespace NICE;

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

SpectralCluster::~SpectralCluster()
{
}

void SpectralCluster::getSimilarityMatrix ( const VVector & features, 
					    // refactor-nice.pl: check this substitution
					    // old: Matrix & laplacian,
					    NICE::Matrix & laplacian,
					    double alpha )
{
    // refactor-nice.pl: check this substitution
    // old: Matrix distances ( laplacian );
    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++ )
    {
	// refactor-nice.pl: check this substitution
	// old: const Vector & xi = features[i];
	const NICE::Vector & xi = features[i];
	for ( int j = i ; j < (int)features.size() ; j++ )
	{
	    // refactor-nice.pl: check this substitution
	    // old: const Vector & xj = features[j];
	    const NICE::Vector & xj = features[j];
	    // double sim = xi * xj;
	    double dist = distance.calculate ( xi, xj );
	    // refactor-nice.pl: check this substitution
	    // old: distances[i][j] = dist;
	    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++ )
	{
	    // refactor-nice.pl: check this substitution
	    // old: double d = ( mean - distances[i][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++ )
	{
	    // refactor-nice.pl: check this substitution
	    // old: double sim = exp(- distances[i][j] * norm );
	    double sim = exp(-distances(i, j) * norm );
	    
	    // refactor-nice.pl: check this substitution
	    // old: laplacian[i][j] = - sim;
	    laplacian(i, j) = - sim;
	    // refactor-nice.pl: check this substitution
	    // old: laplacian[j][i] = - sim;
	    laplacian(j, i) = - sim;
	}
}

void SpectralCluster::computeLaplacian ( const VVector & features,
					 // refactor-nice.pl: check this substitution
					 // old: Matrix & laplacian, int method, double alpha )
					 NICE::Matrix & laplacian, int method, double alpha )
{
	//unused:
    //int n = (int)features.size();
    laplacian.set(0.0);
    getSimilarityMatrix(features, laplacian, alpha);
    
    // refactor-nice.pl: check this substitution
    // old: Vector d ( laplacian.rows() );
    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 VVector & features,
				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);
    }

    // refactor-nice.pl: check this substitution
    // old: const Vector & x = features[0];
    const NICE::Vector & x = features[0];
    int dimension = x.size();

    // refactor-nice.pl: check this substitution
    // old: Matrix laplacian ( features.size(), features.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;
    NICE::eigenvectorvalues ( laplacian, eigvect, eigvals );

    // refactor-nice.pl: check this substitution
    // old: Matrix eigvals_sorted ( eigvals.size(), 1 );
    std::map<double, int> eigvals_sorted;
	for ( int i = 0 ; i < (int)eigvals.size(); i++ ) 
    {
	// refactor-nice.pl: check this substitution
	// old: eigvals_sorted[i][0] = i;
	eigvals_sorted.insert ( make_pair( eigvals[i], i ) );
    }
    // refactor: eigvals_sorted = eigvals_sorted || eigvals;
    // refactor: eigvals_sorted.Sort(1);

    VVector spectral_features;

	for ( int i = 0 ; i < (int)eigvect.rows() ; i++ )
    {
	// refactor-nice.pl: check this substitution
	// old: Vector eigvec_k ( noClasses );
	NICE::Vector eigvec_k ( noClasses );
	map<double, int>::const_iterator k = eigvals_sorted.begin();
	for ( int j = 0 ; j < noClasses ; j++ )
	{
	    int eigval_index = k->second;
	    // refactor-nice.pl: check this substitution
	    // old: eigvec_k[j] = eigvect[i][ (int)eigvals_sorted[j][0] ] ;
	    eigvec_k[j] = eigvect(i, eigval_index ) ;
	    k++;
	}

	spectral_features.push_back ( eigvec_k );
    }

    kmeans.cluster ( spectral_features, prototypes, weights, assignment ); 

    // recompute prototypes

    for ( int k = 0 ; k < noClasses ; 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];
	
	// refactor-nice.pl: check this substitution
	// old: Vector & p = prototypes[k];
	NICE::Vector & p = prototypes[k];
	// refactor-nice.pl: check this substitution
	// old: const Vector & x = *i;
	const NICE::Vector & x = *i;
	p += x;
	weights[k]++;
    }

    for ( int k = 0 ; k < noClasses ; k++ )
    {
	// refactor-nice.pl: check this substitution
	// old: Vector & p = prototypes[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();
    }
}