/*
 * sparse_matrix.cpp
 *
 *  Created on: Nov 18, 2014
 *      Author: qimao
 */

#include "sparse_matrix.h"
#include <utility> // std::swap
#include <armadillo>

using namespace arma;

SeqSparseMatrix::SeqSparseMatrix(Aligner *aligner_, seqan::StringSet<seqan::DnaString> &seqs_)
		:aligner(aligner_),seqs(seqs_)
{
//	similarity_matrix = new dok_matrix();
	num_of_seqs = length(seqs);
	count = 0;
}

SeqSparseMatrix::~SeqSparseMatrix(){}

/*
 * NOTE: in our case, the similarity_matrix is symmetric, so we store the lower diagonal matrix
 * To achieve this, we require i <= j, and the diagonal element is one
 */
double SeqSparseMatrix::get(int i, int j)
{
	if (i > j) std::swap(i,j);

	std::pair<int,int> one_pair(i,j);
	if(similarity_matrix.find(one_pair) == similarity_matrix.end()) // not computed yet, then compute
	{
		float sim_score = 1.0 - aligner->pairwise_distance(seqs[i],seqs[j]);
		similarity_matrix[one_pair] = sim_score;
	}

	return similarity_matrix[one_pair];
}

/*
 * rows and cols should be the global indexes in terms of seqs
 * submatrix is a submatrix in terms of rows and cols
 * openMP parallelizes this section
 */
void SeqSparseMatrix::get_submatrix(uvec &rows, uvec &cols, fmat &submatrix)
{
	int row_len = rows.n_rows;
	int col_len = cols.n_rows;
	submatrix.resize(row_len,col_len);

	for(int i=0;i<row_len;i++)
	{
        #pragma omp parallel for
		for(int j=0;j<col_len;j++)
		{
			int r = rows[i];
			int c = cols[j];

			if (r > c) std::swap(r,c);

			std::pair<int,int> one_pair(r,c);
			if(similarity_matrix.find(one_pair) == similarity_matrix.end())
			{
				float sim_score = 1.0 - aligner->pairwise_distance(seqs[r],seqs[c]);

				#pragma omp critical
				{
					similarity_matrix[one_pair] = sim_score;
					count++;
					if (count % 10000 == 0) { printf("."); fflush(stdout);}
					if (count % 1000000 == 0) printf(" %d \n",count);
				}
				submatrix(i,j) = sim_score;

			}else{
				submatrix(i,j) = similarity_matrix[one_pair];
			}
		}
	}
}

void SeqSparseMatrix::get_submatrix(uword idx, uvec &cols, fvec &submatrix)
{
	int col_len = cols.n_rows;
	submatrix.resize(col_len);

	#pragma omp parallel for
	for(int j=0;j<col_len;j++)
	{
		int r = idx;
		int c = cols[j];

		if (r > c) std::swap(r,c);

		std::pair<int,int> one_pair(r,c);
		if(similarity_matrix.find(one_pair) == similarity_matrix.end())
		{
			float sim_score = 1.0 - aligner->pairwise_distance(seqs[r],seqs[c]);

			#pragma omp critical
			{
				similarity_matrix[one_pair] = sim_score;
				count++;

				if (count % 10000 == 0) { printf("."); fflush(stdout);}
				if (count % 1000000 == 0) printf(" %d \n",count);
			}
			submatrix(j) = sim_score;

		}else{
			submatrix(j) = similarity_matrix[one_pair];
		}
	}
}

