/*
 * test_1.cpp
 *
 *  Created on: Nov 20, 2014
 *      Author: qimao
 */

#include <tclap/CmdLine.h>

#include <seqan/sequence.h>
#include <seqan/stream.h>
#include <seqan/seq_io.h>

#include <iostream>
#include <fstream>
#include <armadillo>
#include <cstdlib>
#include <cstdio>

#include "aligner.h"
#include "flat_cluster.h"
#include "sparse_matrix.h"
#include "hier_cluster.h"
#include "esprit_tree.h"

using namespace arma;

void test_hier_cluster(int argc, char *argv[])
{
	std::string cur_path(argv[0]);
	int last_pos = cur_path.find_last_of("/");
	std::string exe_path = cur_path.substr(0,last_pos);

	Params params(argc,argv);

	std::string fasta_file = params.get_input_filename();
	const char *filename = fasta_file.c_str();
	seqan::SequenceStream seqIO(filename);
	seqan::StringSet<seqan::CharString> ids;
	seqan::StringSet<seqan::DnaString> seqs;

	int res = readAll(ids, seqs, seqIO);
	if (res != 0){
		std::cerr << "ERROR: Could not read records!\n";
		exit(1);
	}

	int distance_opt = params.get_distance_option();
	Aligner * paligner = NULL;
	switch (distance_opt)
	{
	case 0:
		paligner = new KmerAligner(6);
		break;

	case 1:
		paligner = new BandNWAligner(0.1);
		break;

	default:
		std::cout<<"the distance is not defined\n";
		exit(1);
	}

	// apply on the whole data
	int num = length(seqs);
	SeqSparseMatrix smat(paligner,seqs);

	SpectralCluster sc(params);

	ActiveFlatCluster_HalfN afc(params,sc);
	DivisiveHierCluster dhc(afc,params);

	printf("the number of sequence: %d\n",num);

	double t = omp_get_wtime();

	// initialize the number of threads
	int num_proc = params.get_numofthreads();
	omp_set_dynamic(0);     // Explicitly disable dynamic teams
	omp_set_num_threads(num_proc); // Use 4 threads for all consecutive parallel regions

	// perform hybridHC method
	dhc.cluster(smat);

	double t1 = omp_get_wtime() - t;
	printf("\n+++++++AHDC clustering finished! time cost=%f++++++\n",t1);

	// obtain clusters
	std::vector<std::unordered_set<int> > clusters;
	dhc.get_clusters(clusters);

	// apply ESPRIT-Tree method on each clusters
	std::string output_path = params.get_output_filename();
	esprit_tree(clusters,seqs,fasta_file,num_proc,output_path,exe_path);

	t = omp_get_wtime() - t;
	printf ("It took me %f seconds.\n",t);

	int count = smat.get_count();
	printf("The number of pairwise computation: %d\n",count);

	if(output_path != "")
	{
		std::ofstream fout(params.get_output_filename().append("/AHDC.txt"));
		fout<<"the number of seqences: "<<num<<std::endl;
		fout<<"the number of pairwise computation: "<<count<<std::endl;
		fout<<"AHDC time cost: "<<t1<<std::endl;
		fout<<"Total time cost:" << t <<std::endl;
		dhc.print_hier_tree(fout);
		fout.close();
	}else
	{
		std::ostream &fout = std::cout;
		fout<<"the number of seqences: "<<num<<std::endl;
		fout<<"the number of pairwise computation: "<<count<<std::endl;
		fout<<"time cost: "<<t1<<std::endl;
		fout<<"Total time cost:" << t <<std::endl;
		dhc.print_hier_tree(fout);
	}

	delete paligner;
}

int main(int argc, char** argv)
{
	srand(0);
	arma_rng::set_seed(0);
	test_hier_cluster(argc,argv);
//	test_esprit_tree();
}


