#include <stdio.h>
#include "global.h"

REAL	Global::gap_o=10;
REAL	Global::gap_e=0.5;
int	Global::Kmer_Len=5;
REAL *Global::krate=NULL;
int Global::num_krate=1;

REAL Global::krate_step=1.0;
Kmer *Global::kmer=NULL;
DistCache *Global::cache=NULL;
		
Needle *Global::needle=NULL;
Tree ** Global::childlist=NULL;
KDist *Global::kdistlist=NULL;
char *Global::buf1=NULL;
char *Global::buf2=NULL;	
char *Global::ProbSeqOut=NULL;
REAL Global::level_min=0.01;
REAL Global::level_inc=1.25;
REAL Global::level_max=0.15;
REAL Global::level_step=0.01;
int Global::Dist_Style=DS_Ave;
int Global::Needle_Diag=1;
REAL Global::DiagRate=0.1;
int Global::CacheSize=1000;
REAL Global::ErrorRate=1e-8;
REAL Global::currentlevel=0.0;
bool Global::UseApprox=false;
bool Global::ShowAlign=false;

void Global::SetParam(REAL g_o,	REAL g_e, int K_Len, int K_Br)
{
	gap_o=g_o;
	gap_e=g_e;
	Kmer_Len=K_Len;
}

void Global::Init(int SeqNum)
{
	kmer=new Kmer(Kmer_Len);
	needle=new Needle(gap_o,gap_e);
	childlist=(Tree **)Malloc(SeqNum*sizeof(Tree *));
	kdistlist= (KDist *)Malloc(SeqNum*sizeof(KDist));
	buf1=(char *)Malloc(Max_StrLen*sizeof(char));
	buf2=(char *)Malloc(Max_StrLen*sizeof(char));
	cache=new DistCache(CacheSize);
	
}
void Global::Finish()
{
	delete kmer;
	delete needle;
	delete cache;
	free(krate);
	free(childlist);	
	free(kdistlist);
	free(buf1);
	free(buf2);
}

void Global::LoadKList(char *filename)
{
	FILE *fp;
	if (filename==NULL)
	{
		num_krate=1;
		krate=(REAL *)Malloc(sizeof(REAL));
		krate[0]=6.0;
		krate_step=1.0;
		return;		
	}
	if ((fp=fopen(filename,"r"))==NULL)
	{
		fprintf(stderr,"Cannot open Kmer configure file.\n");
		num_krate=1;
		krate=(REAL *)Malloc(sizeof(REAL));
		krate[0]=6.0;
		krate_step=1.0;
	}
	else
	{
		fscanf(fp,"%d",&num_krate);
		fscanf(fp,"%f",&krate_step);
		krate=(REAL *)Malloc(num_krate*sizeof(REAL));
		for (int i=0;i<num_krate;i++)
			fscanf(fp,"%f",krate+i);
		fclose(fp);
	}
}

int CompDist(const void * a, const void * b)
{
	KDist *kd1,*kd2;
	kd1=(KDist *)a;
	kd2=(KDist *)b;
	if (kd1->dist > kd2->dist) return 1;
	if (kd1->dist < kd2->dist) return -1;
	return 0;
}

REAL KdistBound(REAL dist)
{
			int coefidx = (int)(dist/Global::krate_step);
			if (coefidx <0) coefidx=0;
			if (coefidx >=Global::num_krate) coefidx=Global::num_krate-1;
			return Global::krate[coefidx] * dist;
}

REAL calcDistance(char *seq1, char *seq2) {


  unsigned int i, j, k;
  REAL residuecount, distance;
	
  residuecount = distance = 0.0;

	int gap_alert1 = 0;
	int gap_alert2 = 0;

	int start_seq1 = 0;
  int start_seq2 = 0;
  int end_seq1 = strlen(seq1);
  int end_seq2 = strlen(seq2);

	k = 0;
  while(seq1[k] == '-') k++;
  start_seq1 =  k;

  k = 0;
  while(seq2[k] == '-') k++;
  start_seq2 = k;

	k = end_seq1-1;
  while(seq1[k] == '-') k--;
  end_seq1 = k+1; 
  	
	k = end_seq2-1;
  while(seq2[k] == '-') k--;
  end_seq2 = k+1; 
	
  int start_seq = start_seq1;
  if(start_seq < start_seq2)
     {
        start_seq = start_seq2;
     }

  int end_seq = end_seq1;
  if(end_seq > end_seq2)
     {
         end_seq = end_seq2;
     }
	gap_alert1 = 0;
	gap_alert2 = 0;

	for( k=start_seq; k <end_seq; k++) {

	   if (seq1[k] == '-' && seq2[k] == '-')
		{
			continue;
		}

  	if(seq1[k] == 'N' && seq2[k] == 'N')
		{
			gap_alert1 = 0;
			gap_alert2 = 0;
			continue;
		}


  	if( seq1[k] == '.' &&seq2[k] == '.')
		{
			gap_alert1 = 0;
			gap_alert2 = 0;
			continue;
		}


/* multiple consecutive gaps  are counted as 1 mismatch */
/*  but residue is shortened. */
/* you have to check both sequences... */


    if ( seq1[k] == '-' )
		{
			if(gap_alert1 == 0)
			{
        	residuecount+=1.0;
         	distance+=1.0;
			}
			gap_alert1 = 1;
      gap_alert2 = 0;
			continue;
		}

    if ( seq2[k] == '-' )
		{
			if(gap_alert2 == 0)
			{
        	 	residuecount+=1.0;
           	distance+=1.0;
			}
			gap_alert1 = 0;
      gap_alert2 = 1;
			continue;
		}

     /* ok no gaps , check for mismatch */

    if ( seq1[k] !=seq2[k])
		{
    	 distance += 1.0;
       residuecount+=1.0;
			 gap_alert1 = 0;
			 gap_alert2 = 0;
	 		continue;
	 	}

   if ( seq1[k] ==seq2[k])
		{
      residuecount+= 1.0;
			gap_alert1 = 0;
			gap_alert2 = 0;
			continue;	
		}
  } /* end of loop on k */ 

  if (residuecount > 0) {
		distance = distance / residuecount;
  }
  else {
		distance = 1.0;
  }
  return distance;  
}

bool QuickDist(int id1,int id2, int query,int newid, int Sq1, int Sq2, REAL &dist)
{
	if (Global::cache->Get(query,newid,dist)) return true;
	if (id1<0 || id2<0) return false;
	REAL dist1,dist2;
	if (!Global::cache->Get(id1,query,dist1)) return false;
	if (!Global::cache->Get(id2,query,dist2)) return false;
	if (Global::Dist_Style == DS_Ave)
	{
		dist=(Sq1*dist1+Sq2*dist2)/(Sq1+Sq2);
	}
	else
	{
		dist=max(dist1,dist2);
	}
	Global::cache->Store(query,newid,dist);
	return true;
}
