/*
 * BandNWAlign.cpp
 *
 *  Created on: Nov 16, 2014
 *      Author: qimao
 */

#include "aligner.h"

#include <unordered_map>
#include <iostream>

using namespace seqan;

typedef String<Dna> TSequence;
typedef Align<TSequence, ArrayGaps> TAlign;
typedef Row<TAlign>::Type TRow;
typedef Iterator<TRow>::Type TRowIterator;

using namespace std;

int const BandNWAligner::SCORE_MATRIX[TAB_SIZE]={
		 5,-4, -4, -4,
		-4, 5, -4, -4,
		-4, -4, 5, -4,
		-4, -4, -4, 5
	};

//int const BandNWAligner::SCORE_MATRIX[TAB_SIZE]={
//   10, -3, -1, -4,
//   -3, 9, -5, 0,
//   -1, -5, 7, -3,
//   -4, 0, -3, 8
//};

BandNWAligner::BandNWAligner(double band_ratio) {
    this->band_ratio = band_ratio;

    // initialize scoring matrix
    this->scoringSchema = new TScoringSchema(gap_extend_score,gap_open_score);
	for(unsigned i=0;i<ValueSize<Dna>::VALUE; ++i){
		for(unsigned j=0;j<ValueSize<Dna>::VALUE;++j){
			setScore(*scoringSchema,Dna(i),Dna(j), SCORE_MATRIX[i*ValueSize<Dna>::VALUE+j]);
		}
	}
}

double BandNWAligner::pairwise_distance(DnaString &seqH, DnaString &seqV){
	AlignConfig<> alignConfig;

    Align<DnaString> align;
    resize(rows(align), 2);
    assignSource(row(align, 0), seqH);
    assignSource(row(align, 1), seqV);

    // here the diagonal is f(x)=-x, and ensure the last element is inside the band
    int len1 = length(seqH);
    int len2 = length(seqV);
    int lDiag = len2;
    int uDiag = len1;
    if(len1 > len2){
    	int tmp = floor(len2 *band_ratio);
    	uDiag = (len1-len2) + tmp;
    	lDiag = -tmp;
    }else{
    	uDiag = floor(len1 * band_ratio);
    	lDiag =  - (len2-len1 + uDiag);
    }

    //TValue result =
	globalAlignment(align, *scoringSchema, alignConfig, lDiag, uDiag);
//	std::cout<<align;

    // calculate the distance by dropping the gaps in the front and end pairs
    TRow &row1 = row(align,0);
    TRow &row2 = row(align,1);

//  old-inaccurate denomination
//    int left_max = std::max(toViewPosition(row1, 0), toViewPosition(row2, 0));
//    int right_min = std::min(toViewPosition(row1, length(source(row1))-1), toViewPosition(row2, length(source(row2))-1));
//    int valid_len = right_min - left_max + 1;

    std::unordered_map<int,int> set1;
    for(unsigned int i=0;i<length(source(row1)); ++i){
    	set1[toViewPosition(row1, i)] = i;
    }

    double count = 0.0;
    for(unsigned int i=0;i<length(source(row2)); ++i){
    	int idx = toViewPosition(row2, i);
    	if (set1.find(idx) != set1.end() && seqV[i] == seqH[set1[idx]])
    		count++;
    }

    int valid_len = min(len1,len2);
    double dist = 1 - count / valid_len;
    return dist;
}


double KmerAligner::pairwise_distance(DnaString &seqH, DnaString &seqV)
{
	string str_seq1 = "";
	for( Iterator<DnaString>::Type it=seqan::begin(seqH); it !=seqan::end(seqH); ++it)
	{
		char c = *it;
		str_seq1.push_back(c);
	}

	string str_seq2 = "";
	for( Iterator<DnaString>::Type it=seqan::begin(seqV); it !=seqan::end(seqV); ++it)
	{
		char c = *it;
		str_seq2.push_back(c);
	}

	unordered_map<string, int> kmer2count;

	for (unsigned int i=0;i<str_seq1.length()-kmer_window_size+1;++i){
		string kmer = str_seq1.substr(i,kmer_window_size);
		if (kmer2count.find(kmer) == kmer2count.end()){ // does not exist in kmer2count
			kmer2count[kmer] = 1;
		}else{
			kmer2count[kmer] += 1;
		}
	}

/* test map result */
//	unordered_map<string, int>::iterator it = kmer2count.begin();
//	while(it != kmer2count.end()){
//		cout<<"(" << it->first <<","<<it->second<<")"<<endl;
//		it++;
//	}

	int total = 0;
	for(unsigned int i=0;i<str_seq2.length();++i){
		string kmer = str_seq2.substr(i,kmer_window_size);
		if(kmer2count.find(kmer) != kmer2count.end()){
			if(kmer2count[kmer] > 0){
				total += 1;
				--kmer2count[kmer];
			}
		}
	}

	double dist = 1.0 - 1.0 * total / ( min(str_seq1.length(), str_seq2.length()) - kmer_window_size + 1 );

	return dist;
}

