/*
 * flatclusteralg.cpp
 *
 *  Created on: Nov 19, 2014
 *      Author: qimao
 */

#include "flat_cluster.h"

#include <vector>
#include <cfloat>
#include <cmath>
#include <cstdlib>

float KMeans::cluster_rand_init(fmat & W, fmat &centroids_, uvec &clusters_)
{
	int num_sample = W.n_rows;
	int dim = W.n_cols;
	centroids_.resize(k,dim);
	clusters_.resize(num_sample);

	// randomly initializing centroids
	for(int c=0;c<k;c++) centroids_.row(c) = W.row( rand() % num_sample );

	int iter = 0;
	float old_mse = FLT_MAX;
	while (iter < maxiter)
	{
		//assignment
		fvec sum_distsq = zeros<fvec>(k);
		uvec count = zeros<uvec>(k);
		for(int i = 0;i<num_sample;++i)
		{
			int min_idx=-1;
			float min_dist=FLT_MAX;
			frowvec vec_i = W.row(i);
			for(int c=0;c<k;++c)
			{
				frowvec vec_c = centroids_.row(c);
				float dist = get_distance(vec_i,vec_c);
				if (c==0 || min_dist > dist)
				{
					min_dist = dist;
					min_idx = c;
				}
			}
			clusters_(i) = min_idx;
			sum_distsq(min_idx) += min_dist * min_dist;
			count(min_idx)++;
		}

		// mean square error
		float mse = 0.0;
		for(int c=0;c<k;++c) mse += sum_distsq(c);

//		std::cout<<"iter="<<iter<<", mse="<<mse<<std::endl;

		if (iter > 1 && (old_mse - mse)/old_mse < this->relative_mse)
		{
			old_mse = mse;
			break;
		}

		// update the centroids
		centroids_.zeros(k,dim);

		for(int i =0;i<num_sample;++i)
			centroids_.row(clusters_(i)) += W.row(i) / count(clusters_(i));

		iter++;
		old_mse = mse;
	}

//	std::cout<<"the current mse="<<old_mse<<std::endl;
	return old_mse;
}

void KMeans::cluster(fmat &W)
{
	int num_sample = W.n_rows;
	int dim = W.n_cols;
	centroids.resize(k,dim);
	clusters.resize(num_sample);

	float best_mse = FLT_MAX;
	for(int i = 0;i<nrepeats;i++)
	{
//		std::cout<<"the "<<i<<"th repeat"<<std::endl;

		fmat centroids_;
		uvec clusters_;
		float mse = cluster_rand_init(W, centroids_,clusters_);

		if(i == 0 || best_mse > mse)
		{
			best_mse = mse;
			for(int j=0;j<num_sample;++j) clusters(j) = clusters_(j);
			for(int c=0;c<k;++c) centroids.row(c) = centroids_.row(c);
		}
	}
//	std::cout<<"the best mse is "<<best_mse<<std::endl;
}

void SpectralCluster::cluster(fmat & W)
{
	fmat D = diagmat( sum(W,1) );
	fmat L = D - W;
	fvec eigval;
	fmat eigvec;
	eig_sym(eigval, eigvec, L);
	uvec index = sort_index(eigval,"ascend");

	int num_sample = W.n_rows;
	fmat M(num_sample,k);
	for(int col=0;col<k;++col)
		M.col(col) = eigvec.col(index(col));

	KMeans kmeans(params);
	kmeans.cluster(M);

	uvec p = kmeans.get_clusters();
	clusters.resize(num_sample);
	for(int i=0;i<num_sample;++i) clusters(i) = p(i);
}


void ActiveFlatCluster::landmark_selection(SeqSparseMatrix &smat, uvec &indices, uvec &sub_indices)
{
	int n = indices.n_rows;
	unsigned int s = floor( log2 (n * 1.0) ); // fixing the name of samples used

	sub_indices.resize(s);
	std::unordered_set<int> landmark_set; // maintain a set of unique indices
	int i = 0;
	while (landmark_set.size() < s) // terminate when landmark_set = s
	{
		int idx = indices( rand() % n );
		if (landmark_set.find(idx) == landmark_set.end())
		{
			landmark_set.insert(idx);
			sub_indices(i)=idx;
			i++;
		}
	}

	// compute the approximate diameter
	fmat W;
	smat.get_submatrix(sub_indices,sub_indices,W);
	this->approx_diameter_dist = 1.0 - W.min();

	printf("random landmark selection, diameter=%f\n",this->approx_diameter_dist);
}

void ActiveFlatCluster::active_cluster(SeqSparseMatrix &smat, uvec & indices)
{
	// reset clusters
	clusters.clear(); // remember to clear all clusters for the following calls
	for(int c=0;c<k;++c)
	{
		std::unordered_set<uword> one_set;
		clusters.push_back(one_set);
	}

	// select subset as landmarks
	uvec sub_indices;
	landmark_selection(smat,indices,sub_indices);

	// perform flat clustering on the landmark set
	fmat W;
	smat.get_submatrix(sub_indices,sub_indices,W);
	base_alg.cluster(W);
	uvec base_clusters = base_alg.get_clusters();

	// assign the landmarks
	fvec count = zeros<fvec>(k);
	std::unordered_set<uword> landmarks;
	for(unsigned int i = 0;i<sub_indices.n_rows;++i)
	{
		int cluster_idx = base_clusters(i);
		int global_idx = sub_indices(i);
		clusters[cluster_idx].insert(global_idx);
		landmarks.insert(global_idx);
		count(cluster_idx)++;
	}

	// assign all the rest of data to clusters
#pragma omp parallel for
	for(unsigned int i =0;i<indices.n_rows;++i)
	{
		int cur_idx = indices(i);
		if(landmarks.find(cur_idx) == landmarks.end()) //not find in the landmark set
		{
			fvec sum_dist = zeros<fvec>(k);
			for(unsigned int j=0;j<base_clusters.n_rows;++j)
			{
				float dist = 1.0 - smat.get(cur_idx,sub_indices(j)); // time consuming part, similarity to distance
				sum_dist(base_clusters(j)) += dist;
			}

			float ave_min = FLT_MAX;
			int min_idx = -1;
			for(int c=0;c<k;++c)
			{
				float tmp = sum_dist(c)/count(c);   // assigned to cluster c with the smallest average distance
				if(c==0 || ave_min > tmp)
				{
					ave_min = tmp;
					min_idx = c;
				}
			}

		#pragma omp critical
			clusters[min_idx].insert(cur_idx);
		}
	}
}

void ActiveFlatCluster_HalfN::landmark_selection(SeqSparseMatrix &smat, uvec & indices, uvec & sub_indices)
{
	int n = indices.n_rows;
	int s = floor( log2 (n * 1.0) ); // fixing the name of samples used
	int q = n / 2;

	sub_indices.resize(s); // keep the middle sampled
	std::unordered_set<int> landmark_set;

	// initially sample one as the start
	int idx = indices( rand() % n );
	sub_indices(0) = idx;
	landmark_set.insert(idx);

	// maintain the maximal similarity of each sequence to the selected set of landmarks
	fvec max_sim;
	smat.get_submatrix(idx,indices,max_sim);

	for(int i = 1; i<s; ++i)
	{
		uvec order_vec = sort_index(max_sim,"ascend");
		idx = indices( order_vec( rand() % (n-q-1) ));
		while( landmark_set.find(idx) != landmark_set.end() )
			idx = indices( order_vec( rand() % (n-q-1) ));
		sub_indices(i) = idx;
		landmark_set.insert(idx);

		// update max_sim for similarity between one seq to a set of selected landmarks
		fvec idx_sim;
		smat.get_submatrix(idx,indices,idx_sim);
	#pragma omp parallel for
		for(int j=0;j<n;++j)
			max_sim(j) = std::max(max_sim(j), idx_sim(j));
	}

// compute the approximate diameter
	fmat W;
	smat.get_submatrix(sub_indices,sub_indices,W);
	this->approx_diameter_dist = 1.0 - W.min();

//	// find the second farthest points
//	// take the advantage of the previous computed pairwise distance
//
//	// use 2 sampled instances for estimate the diameter
//	int tr_num = 3;
//	fvec min_sim(tr_num);
//#pragma omp parallel for
//	for(int i=0; i< tr_num; ++i)
//	{
//		int start_idx = sub_indices(i);
//
//		fvec start_idx_sim;
//		smat.get_submatrix(start_idx,indices,start_idx_sim);
//		uword min_idx;
//		start_idx_sim.min(min_idx);
//		int second_idx = indices(min_idx);
//
//		fvec second_idx_sim;
//		smat.get_submatrix(second_idx,indices,second_idx_sim);
//		min_sim(i) = second_idx_sim.min(min_idx);
//	}
//
//	// compute the approximate diameter
//	this->approx_diameter_dist = 1.0 - min_sim.min();
	printf("half N landmark selection, diameter=%f\n",this->approx_diameter_dist);
}

