#include <stdlib.h>
#include <string.h>
#include <ctype.h>
#include "needle.h"
#include "util.h"
#include <iostream>
#include <fstream>
using namespace std;

#define MinFloat -1e100
                                 // A    G   C   T
static REAL simularity[] = {10, -1, -3, -4, //A
                            -1,  7, -5, -3, //G
                            -3, -5,  9,  0, //C
                            -4, -3,  0,  8};//T    

DIRS VTracks::GetMaxInc(REAL ia, REAL ib, REAL ic, REAL &newval)
{
	DIRS dirs=DIR_UP;
	newval=vals[DIR_UP]+ic;
	if (vals[DIR_LEFT] + ib > newval)
		{
			newval= vals[DIR_LEFT] + ib;
			dirs=DIR_LEFT;
		}	
	if (vals[DIR_DIAG] + ia > newval)
		{
			newval= vals[DIR_DIAG] + ia;
			dirs=DIR_DIAG;
		}	
	return dirs;
}

REAL Needle::Score(ProbChar &p1, ProbChar &p2)
{
	REAL ret=0.0;
	for (int i=1;i<=Num_RNAType;i++)
	{
		ret+=p1[i]*(9*p2[i]-4*(1-p2[0]-p2[Num_RNAType+2]));
	}
	ret-=4*p1[Num_RNAType+1]*(1-p2[0]-p2[Num_RNAType+1]-p2[Num_RNAType+2]);
	return ret;
}

REAL Needle::Score(int row, ProbChar &p2)
{
	if (row>4|| row <=-1) return 0;
	if (row==-1) cerr << "here "<<endl;
	if (row==4) return -4*(1-p2[0]-p2[Num_RNAType+1]-p2[Num_RNAType+2]);
	return 9*p2[row+1]-4*(1-p2[0]-p2[Num_RNAType+2]);
}

REAL Needle::Score(int row, int col)
{
    if (row==4 && col==4) return 0;
    if (row==col)
    	return 5;
   	return -4;
    //return simularity[row * 4 + col];
}


int* Needle::StrToCode(char* str)
{
	int len=strlen(str);
	int *l=(int *)Malloc(strlen(str)*sizeof(int));
	
	for (int i = 0; i < len; i++)
	{
		switch(toupper(str[i]))
		{
		case 'A':
				l[i]=0;
				break;
		case 'G':
				l[i]=1;
				break;
		case 'C':
				l[i]=2;
				break;
		case 'T':
				l[i]=3;
				break;
		case 'N':
			  l[i]=4;
			  break;		
		default:
				l[i]=-1;
				break;
		}
	}
	return l;
}

REAL Needle::CalculateMatrix(int* source, int* dest, int slen, int dlen,DIRTracks *pr)
{
  int i,j,t;
  REAL inc,kdist;
 // FILE *fp;
  
	VTracks *ar = new VTracks[dlen+1];
  VTracks lefttrack,thistrack;
	DIRS dirs;
 
 	//fp=fopen("nul","w");
 	
 	ar[0].SetVal(DIR_DIAG,0);
 	ar[0].SetVal(DIR_LEFT,MinFloat);
 	ar[0].SetVal(DIR_UP,MinFloat);
  
  for (i = 1; i <= dlen; i++)
  {
  	ar[i].SetVal(DIR_DIAG,MinFloat);
  	ar[i].SetVal(DIR_UP,MinFloat);
    ar[i].SetVal(DIR_LEFT,gap_len*i);
  }
        
        
  for (i = 1; i <= slen; i++)
  {
  	lefttrack.SetVal(DIR_DIAG,MinFloat);
  	lefttrack.SetVal(DIR_UP,gap_len*i);
  	lefttrack.SetVal(DIR_LEFT,MinFloat);
  	
    for (j = 1; j <= dlen; j++)
      {  
  			dirs=ar[j-1].GetMaxInc(0,0,0,kdist);
			kdist+=Score(source[i-1],dest[j-1]);                 
			pr[i*(dlen+1)+j].SetVal(DIR_DIAG,dirs);
			thistrack.SetVal(DIR_DIAG,kdist);
			if (i==slen)
				dirs=lefttrack.GetMaxInc(gap_len,gap_len,gap_len,kdist);
			else
				dirs=lefttrack.GetMaxInc(gap_open,gap_len,gap_open,kdist);
			pr[i*(dlen+1)+j].SetVal(DIR_LEFT,dirs);
			thistrack.SetVal(DIR_LEFT,kdist);
  			if (j==dlen)
				dirs=ar[j].GetMaxInc(gap_len,gap_len,gap_len,kdist);
			else
				dirs=ar[j].GetMaxInc(gap_open,gap_open,gap_len,kdist);
			pr[i*(dlen+1)+j].SetVal(DIR_UP,dirs);
			thistrack.SetVal(DIR_UP,kdist);
			ar[j-1]=lefttrack;
			lefttrack=thistrack;
      }
    ar[dlen]=thistrack;  

//	  for (t=1;t<=dlen;t++)
 //	 {
 // 		fprintf(fp,"(%5.1f%5.1f%5.1f)",ar[t].GetVal(DIR_DIAG),ar[t].GetVal(DIR_LEFT),ar[t].GetVal(DIR_UP)); 
 // 	}      
//	  for (t=1;t<=dlen;t++)
//	 	{
//  		fprintf(fp,"(%5d%5d%5d)",pr[i][t].GetVal(DIR_DIAG),pr[i][t].GetVal(DIR_LEFT),pr[i][t].GetVal(DIR_UP)); 
 // 	}      
//  	fprintf(fp,"\n");
  }
  
  dirs=ar[dlen].GetMaxInc(0,0,0,kdist);
  pr[slen*(dlen+1)+dlen].SetEntry(dirs);
//	fprintf(fp,"%d %f\n",dirs,kdist);
//  fclose(fp);

  delete [] ar;
  return kdist;
}
  
REAL Needle::CalculateMatrix(ProbString &source, ProbString &dest, int slen, int dlen, DIRTracks *pr)
{
  int i,j,t;
  REAL inc,kdist;

	VTracks *ar = new VTracks[dlen+1];
	VTracks lefttrack,thistrack;
	DIRS dirs;
	 
	 
 	ar[0].SetVal(DIR_DIAG,0);
 	ar[0].SetVal(DIR_LEFT,MinFloat);
 	ar[0].SetVal(DIR_UP,MinFloat);
 	
 	REAL lastv;
 	lastv=0;//gap_open*(1-dest[0][0]);
	for (i = 1; i <= dlen; i++)
	{
		ar[i].SetVal(DIR_DIAG,MinFloat);
		ar[i].SetVal(DIR_UP,MinFloat);
		ar[i].SetVal(DIR_LEFT,lastv);
		if (i<dlen)
			lastv+=gap_len*(1-dest[i][0]);
	}
        
	lastv=gap_len*(1-source[0][0]);
        
	for (i = 1; i <= slen; i++)
	{
		lefttrack.SetVal(DIR_DIAG,MinFloat);
		lefttrack.SetVal(DIR_UP,lastv);
		lefttrack.SetVal(DIR_LEFT,MinFloat);
		if (i<slen)
			lastv+=gap_len*(1-source[i][0]);
  	
		for (j = 1; j <= dlen; j++)
		{  
			dirs=ar[j-1].GetMaxInc(0,0,0,kdist);
			kdist+=Score(source[i-1],dest[j-1]);
			pr[i*(dlen+1)+j].SetVal(DIR_DIAG,dirs);
			thistrack.SetVal(DIR_DIAG,kdist);
			if (i==slen)
				dirs=lefttrack.GetMaxInc(gap_len*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),kdist);
			else
				dirs=lefttrack.GetMaxInc(gap_open*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),gap_open*(1-dest[j-1][0]),kdist);
			pr[i*(dlen+1)+j].SetVal(DIR_LEFT,dirs);
			thistrack.SetVal(DIR_LEFT,kdist);
			if (j==dlen)
				dirs=ar[j].GetMaxInc(gap_len*(1-source[i-1][0]),gap_len*(1-source[i-1][0]),gap_len*(1-source[i-1][0]),kdist);
			else
				dirs=ar[j].GetMaxInc(gap_open*(1-source[i-1][0]),gap_open*(1-source[i-1][0]),gap_len*(1-source[i-1][0]),kdist);
			pr[i*(dlen+1)+j].SetVal(DIR_UP,dirs);
			thistrack.SetVal(DIR_UP,kdist);
			ar[j-1]=lefttrack;
			lefttrack=thistrack;
		}
		ar[dlen]=thistrack;  
	}
  
	dirs=ar[dlen].GetMaxInc(0,0,0,kdist);
	pr[slen*(dlen+1)+dlen].SetEntry(dirs);
//	fprintf(fp,"%d %f\n",dirs,kdist);
//  fclose(fp);

	delete [] ar;
	return kdist;
}
  
REAL Needle::CalculateMatrix_Diag(int *source, int *dest, int slen, int dlen, int diag_dn, int diag_up, DIRTracks *pr)
{
	int i,j,t;
	REAL inc,kdist;

	VTracks *ar = new VTracks[dlen+1];
	VTracks lefttrack,thistrack;
	DIRS dirs;
	
	 
	int rdif=max(dlen-slen,0);
	int ldif=max(slen-dlen,0);
	 
 	ar[0].SetVal(DIR_DIAG,0);
 	ar[0].SetVal(DIR_LEFT,MinFloat);
 	ar[0].SetVal(DIR_UP,MinFloat);
 	
 	REAL lastv;
 	lastv=0;
	for (i = 1; i <= dlen; i++)
	{
  		ar[i].SetVal(DIR_DIAG,MinFloat);
  		ar[i].SetVal(DIR_UP,MinFloat);
		if (i<=rdif+diag_up)
		{
			ar[i].SetVal(DIR_LEFT,lastv);
			if (i<dlen)
				lastv+=gap_len;
		}
		else
		{
			ar[i].SetVal(DIR_LEFT,MinFloat);
		}
	}
        
 	lastv=0;//gap_open;
        
	for (i = 1; i <= slen; i++)
	{
		lefttrack.SetVal(DIR_DIAG,MinFloat);
		lefttrack.SetVal(DIR_LEFT,MinFloat);
		if (i<=ldif+diag_dn)
		{
			lefttrack.SetVal(DIR_UP,lastv);
			if (i<slen)
				lastv+=gap_len;
		}
		else
		{
			lefttrack.SetVal(DIR_UP,MinFloat);
		}
		int colstart=max(1,i-ldif-diag_dn);
 		int colend=min(dlen,i+rdif+diag_up);
		for (j = colstart; j <= colend; j++)
		{  
  			dirs=ar[j-1].GetMaxInc(0,0,0,kdist);
			kdist+=Score(source[i-1],dest[j-1]);
			pr[i*(dlen+1)+j].SetVal(DIR_DIAG,dirs);
			thistrack.SetVal(DIR_DIAG,kdist);
			if (i==slen)
				dirs=lefttrack.GetMaxInc(gap_len,gap_len,gap_len,kdist);
			else
				dirs=lefttrack.GetMaxInc(gap_open,gap_len,gap_open,kdist);
			pr[i*(dlen+1)+j].SetVal(DIR_LEFT,dirs);
			thistrack.SetVal(DIR_LEFT,kdist);
			if (j==dlen)
				dirs=ar[j].GetMaxInc(gap_len,gap_len,gap_len,kdist);
			else
				dirs=ar[j].GetMaxInc(gap_open,gap_open,gap_len,kdist);
			pr[i*(dlen+1)+j].SetVal(DIR_UP,dirs);
			thistrack.SetVal(DIR_UP,kdist);
			ar[j-1]=lefttrack;
			lefttrack=thistrack;
		}
		ar[colend]=thistrack;  
	}
  
	dirs=ar[dlen].GetMaxInc(0,0,0,kdist);
	pr[slen*(dlen+1)+dlen].SetEntry(dirs);
//	fprintf(fp,"%d %f\n",dirs,kdist);
//  fclose(fp);

	delete [] ar;
	return kdist;
}

REAL Needle::CalculateMatrix_Diag(ProbString &source, ProbString &dest, int slen, int dlen, int diag_dn, int diag_up, DIRTracks *pr)
{
  int i,j,t;
  REAL inc,kdist;

	VTracks *ar = new VTracks[dlen+1];
	VTracks lefttrack,thistrack;
	DIRS dirs;
	
	 
	int rdif=max(dlen-slen,0);
	int ldif=max(slen-dlen,0);
	 
 	ar[0].SetVal(DIR_DIAG,0);
 	ar[0].SetVal(DIR_LEFT,MinFloat);
 	ar[0].SetVal(DIR_UP,MinFloat);
 	
 	REAL lastv;
 	lastv=0; //gap_open*(1-dest[0][0]);
	for (i = 1; i <= dlen; i++)
	{
  		ar[i].SetVal(DIR_DIAG,MinFloat);
  		ar[i].SetVal(DIR_UP,MinFloat);
		if (i<=rdif+diag_up)
		{
			ar[i].SetVal(DIR_LEFT,lastv);
			if (i<dlen)
				lastv+=gap_len*(1-dest[i][0]);
		}
		else
		{
			ar[i].SetVal(DIR_LEFT,MinFloat);
		}
	}
        
 	lastv=0; //gap_open*(1-source[0][0]);
        
	for (i = 1; i <= slen; i++)
	{
		lefttrack.SetVal(DIR_DIAG,MinFloat);
		lefttrack.SetVal(DIR_LEFT,MinFloat);
		if (i<=ldif+diag_dn)
		{
			lefttrack.SetVal(DIR_UP,lastv);
			if (i<slen)
				lastv+=gap_len*(1-source[i][0]);
		}
		else
		{
			lefttrack.SetVal(DIR_UP,MinFloat);
		}
 		int colstart=max(1,i-ldif-diag_dn);
 		int colend=min(dlen,i+rdif+diag_up);
		for (j = colstart; j <= colend; j++)
		{  
			dirs=ar[j-1].GetMaxInc(0,0,0,kdist);
			kdist+=Score(source[i-1],dest[j-1]);
			pr[i*(dlen+1)+j].SetVal(DIR_DIAG,dirs);
			thistrack.SetVal(DIR_DIAG,kdist);
			if (i==slen)
				dirs=lefttrack.GetMaxInc(gap_len*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),kdist);
			else
				dirs=lefttrack.GetMaxInc(gap_open*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),gap_open*(1-dest[j-1][0]),kdist);
			pr[i*(dlen+1)+j].SetVal(DIR_LEFT,dirs);
			thistrack.SetVal(DIR_LEFT,kdist);
			if (j==dlen)
				dirs=ar[j].GetMaxInc(gap_len*(1-source[i-1][0]),gap_len*(1-source[i-1][0]),gap_len*(1-source[i-1][0]),kdist);
			else
				dirs=ar[j].GetMaxInc(gap_open*(1-source[i-1][0]),gap_open*(1-source[i-1][0]),gap_len*(1-source[i-1][0]),kdist);
			pr[i*(dlen+1)+j].SetVal(DIR_UP,dirs);
			thistrack.SetVal(DIR_UP,kdist);
			ar[j-1]=lefttrack;
			lefttrack=thistrack;
		}
		ar[colend]=thistrack;  
	}
  
	dirs=ar[dlen].GetMaxInc(0,0,0,kdist);
	pr[slen*(dlen+1)+dlen].SetEntry(dirs);
//	fprintf(fp,"%d %f\n",dirs,kdist);
//  fclose(fp);

	delete [] ar;
	return kdist;
}

REAL Needle::CalculateMatrix_Diag(int *source, ProbString &dest, int slen, int dlen, int diag_dn, int diag_up, DIRTracks *pr)
{
	int i,j,t;
	REAL inc,kdist;

	VTracks *ar = new VTracks[dlen+1];
	VTracks lefttrack,thistrack; 
	DIRS dirs;
	
	 
	int rdif=max(dlen-slen,0);
	int ldif=max(slen-dlen,0);
	 
 	ar[0].SetVal(DIR_DIAG,0);
 	ar[0].SetVal(DIR_LEFT,MinFloat);
 	ar[0].SetVal(DIR_UP,MinFloat);
 	
 	REAL lastv;
 	lastv=0; //gap_open*(1-dest[0][0]);
	for (i = 1; i <= dlen; i++)
	{
  		ar[i].SetVal(DIR_DIAG,MinFloat);
  		ar[i].SetVal(DIR_UP,MinFloat);
		if (i<=rdif+diag_up)
		{
			ar[i].SetVal(DIR_LEFT,lastv);
			if (i<dlen)
    		lastv+=gap_len*(1-dest[i][0]);
		}
		else
		{
			ar[i].SetVal(DIR_LEFT,MinFloat);
		}
	}
        
 	lastv=0; //gap_open;
        
	for (i = 1; i <= slen; i++)
	{
		lefttrack.SetVal(DIR_DIAG,MinFloat);
		lefttrack.SetVal(DIR_LEFT,MinFloat);
		if (i<=ldif+diag_dn)
		{
			lefttrack.SetVal(DIR_UP,lastv);
			if (i<slen)
				lastv+=gap_len;
		}
		else
		{
			lefttrack.SetVal(DIR_UP,MinFloat);
 		}
 		int colstart=max(1,i-ldif-diag_dn);
 		int colend=min(dlen,i+rdif+diag_up);
		for (j = colstart; j <= colend; j++)
		{  
			dirs=ar[j-1].GetMaxInc(0,0,0,kdist);
			kdist+=Score(source[i-1],dest[j-1]);
			pr[i*(dlen+1)+j].SetVal(DIR_DIAG,dirs);
			thistrack.SetVal(DIR_DIAG,kdist);
			if (i==slen)
				dirs=lefttrack.GetMaxInc(gap_len*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),kdist);
			else
				dirs=lefttrack.GetMaxInc(gap_open*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),gap_open*(1-dest[j-1][0]),kdist);
			pr[i*(dlen+1)+j].SetVal(DIR_LEFT,dirs);
			thistrack.SetVal(DIR_LEFT,kdist);
			if (j==dlen)
				dirs=ar[j].GetMaxInc(gap_len,gap_len,gap_len,kdist);
			else
				dirs=ar[j].GetMaxInc(gap_open,gap_open,gap_len,kdist);
			pr[i*(dlen+1)+j].SetVal(DIR_UP,dirs);
			thistrack.SetVal(DIR_UP,kdist);
			ar[j-1]=lefttrack;
			lefttrack=thistrack;
		}
		ar[colend]=thistrack;  
	}
  
	dirs=ar[dlen].GetMaxInc(0,0,0,kdist);
	pr[slen*(dlen+1)+dlen].SetEntry(dirs);
//	fprintf(fp,"%d %f\n",dirs,kdist);
//  fclose(fp);

	delete [] ar;
	return kdist;
}

REAL Needle::CalculateMatrix_NoDgap(int *source, ProbString &dest, int slen, int dlen, DIRTracks *pr)
{
	int i,j,t;
	REAL inc,kdist;
	if (dlen < slen)
		cerr <<" Error Template len < string len " <<endl;

	VTracks *ar = new VTracks[dlen+1];
	VTracks lefttrack,thistrack; 
	DIRS dirs;
	
	 
 	ar[0].SetVal(DIR_DIAG,0);
 	ar[0].SetVal(DIR_LEFT,MinFloat);
 	ar[0].SetVal(DIR_UP,MinFloat);
 	
 	REAL lastv;
 	lastv=0; //gap_open*(1-dest[0][0]);
	for (i = 1; i <= dlen; i++)
	{
  		ar[i].SetVal(DIR_DIAG,MinFloat);
  		ar[i].SetVal(DIR_UP,MinFloat);
		if (i==1)
		{
			ar[i].SetVal(DIR_LEFT,lastv);
		}
		else
		{
			ar[i].SetVal(DIR_LEFT,MinFloat);
		}
	}
        
 	lastv=0; //gap_open;
        
	for (i = 1; i <= slen; i++)
	{
		lefttrack.SetVal(DIR_DIAG,MinFloat);
		lefttrack.SetVal(DIR_LEFT,MinFloat);
		lefttrack.SetVal(DIR_UP,lastv);
		for (j = 1; j <= i+(dlen-slen); j++)
		{  
			dirs=ar[j-1].GetMaxInc(0,0,0,kdist);
			kdist+=Score(source[i-1],dest[j-1]);
			pr[i*(dlen+1)+j].SetVal(DIR_DIAG,dirs);
			thistrack.SetVal(DIR_DIAG,kdist);
			if (i==slen)
				dirs=lefttrack.GetMaxInc(gap_len*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),kdist);
			else
				dirs=lefttrack.GetMaxInc(gap_open*(1-dest[j-1][0]),gap_len*(1-dest[j-1][0]),gap_open*(1-dest[j-1][0]),kdist);
			pr[i*(dlen+1)+j].SetVal(DIR_LEFT,dirs);
			thistrack.SetVal(DIR_LEFT,kdist);
			thistrack.SetVal(DIR_UP,MinFloat);
			ar[j-1]=lefttrack;
			lefttrack=thistrack;
		}
		ar[i+(dlen-slen)]=thistrack;  
	}
  
	dirs=ar[dlen].GetMaxInc(0,0,0,kdist);
	pr[slen*(dlen+1)+dlen].SetEntry(dirs);
//	fprintf(fp,"%d %f\n",dirs,kdist);
//  fclose(fp);

	delete [] ar;
	return kdist;
}

void Needle::GetAlignments(DIRTracks *pr, char *sA, char *sB, char *alA,char *alB)
{
  char *buf,*ptrA,*ptrB,*pbuf;
	int i,j;
	DIRS endpt,newpt;
	
  ptrA=alA;
  ptrB=alB;       
         
  i=strlen(sA);
  j=strlen(sB);    
  int dlen=j;

	buf=(char *)Malloc((i+j+1)*sizeof(char));
        
  while (i > 0 && j > 0)
  {
    endpt=pr[i*(dlen+1)+j].GetEntry();
    newpt=pr[i*(dlen+1)+j].GetEntryVal();
    if (endpt==DIR_DIAG)
       {
         *(ptrA++) = sA[i-1];
         *(ptrB++) = sB[j-1];
         i--;j--;                
       }
    else if (endpt==DIR_UP)
      {
        *(ptrA++) = sA[i-1];
        *(ptrB++) = '-';
        i--;
      }
    else if (endpt==DIR_LEFT)
      {
        *(ptrA++) = '-';
        *(ptrB++) = sB[j-1];
        j--;
      }
    else
    {	
    	fprintf(stderr,"Error No Match %d %d\n",i,j);
    	exit(0);
    }  
    pr[i*(dlen+1)+j].SetEntry(newpt);
  }
  while(i > 0)
  {
    *(ptrA++) = sA[i-1];
    *(ptrB++) = '-';
    i--;            
  }
  while(j > 0)
  {
    *(ptrA++) = '-';
    *(ptrB++) = sB[j - 1];
    j--;            
  }
	ptrA--;ptrB--;

	pbuf=buf;
	do
	{
		*(pbuf++)=*(ptrA--);
	}while (ptrA >=alA);
	*(pbuf)=0;
	strcpy(alA,buf);
	  
	pbuf=buf;
	do
	{
		*(pbuf++)=*(ptrB--);
	}while (ptrB >=alB);
	*(pbuf)=0;
	strcpy(alB,buf);
	free (buf);
}

void Needle::GetAlignments(DIRTracks *pr, ProbString &sA, ProbString &sB, ProbString &alA,ProbString &alB)
{
	ProbChar *bufA,*bufB,*ptrA,*ptrB;
	int i,j;
	DIRS endpt,newpt;

	i=sA.Len();
	j=sB.Len();          
	
	int dlen=j;
	
	bufA=new ProbChar[i+j];
	bufB=new ProbChar[i+j];
	
	ptrA=bufA;
	ptrB=bufB;       
       
	while (i > 0 && j > 0)
	{
		endpt=pr[i*(dlen+1)+j].GetEntry();
		newpt=pr[i*(dlen+1)+j].GetEntryVal();
		if (endpt==DIR_DIAG)
		{
			*(ptrA++) = sA[i-1];
			*(ptrB++) = sB[j-1];
			i--;j--;                
		}
		else if (endpt==DIR_UP)
		{
			*(ptrA++) = sA[i-1];
			*(ptrB++) = gapChar;
			i--;
		}
		else if (endpt==DIR_LEFT)
		{
			*(ptrA++) = gapChar;
			*(ptrB++) = sB[j-1];
			j--;
		}
		else
		{	
			cerr <<"Error No Match" << i << " " << j << endl;
			cerr <<sA.ToString()<<endl;
			cerr <<sB.ToString()<<endl;
			exit(0);
		}  
		pr[i*(dlen+1)+j].SetEntry(newpt);
	}
	while(i > 0)
	{
		*(ptrA++) = sA[i-1];
		*(ptrB++) = gapChar;
		i--;            
	}
	while(j > 0)
	{
		*(ptrA++) = gapChar;
		*(ptrB++) = sB[j - 1];
		j--;            
	}
	
	
	alA.Alloc((int) (ptrA-bufA));
	alB.Alloc((int) (ptrB-bufB));
	
	for (i=0;i<alA.Len();i++) alA[i]=*(--ptrA);
	for (i=0;i<alB.Len();i++) alB[i]=*(--ptrB);
	if (ptrA!=bufA || ptrB !=bufB)
	{
		cerr << "Needle align ptr error\n"<<endl;
	}
	//convert gap to end-gap
	int ptr;
	
	i=0;
	while (i <alA.Len() && alA[i].isGap()) alA[i++]=EndgapChar;
	j=alA.Len()-1;
	while (j>=0 && alA[j].isGap()) 	alA[j--]=EndgapChar;
	for (ptr=i+1;ptr<j;ptr++)
	{
		if (alA[ptr].isGap())
		{
			alA[ptr][Num_RNAType+2]=max(alA[ptr-1][Num_RNAType+2],alA[ptr][Num_RNAType+2]);
			alA[ptr][0]=1.0-alA[ptr][Num_RNAType+2];
		}
	}
	for (ptr=j-1;ptr>i;ptr--)
	{
		if (alA[ptr].isGap())
		{
			alA[ptr][Num_RNAType+2]=max(alA[ptr+1][Num_RNAType+2],alA[ptr][Num_RNAType+2]);
			alA[ptr][0]=1.0-alA[ptr][Num_RNAType+2];
		}
	}
	
	i=0;
	while (i <alB.Len() && alB[i].isGap()) alB[i++]=EndgapChar;
	j=alB.Len()-1;
	while (j>=0 && alB[j].isGap()) alB[j--]=EndgapChar;
	for (ptr=i+1;ptr<j;ptr++)
	{
		if (alB[ptr].isGap())
		{
			alB[ptr][Num_RNAType+2]=max(alB[ptr-1][Num_RNAType+2],alB[ptr][Num_RNAType+2]);
			alB[ptr][0]=1.0-alB[ptr][Num_RNAType+2];
		}
	}
	for (ptr=j-1;ptr>i;ptr--)
	{
		if (alB[ptr].isGap())
		{
			alB[ptr][Num_RNAType+2]=max(alB[ptr+1][Num_RNAType+2],alB[ptr][Num_RNAType+2]);
			alB[ptr][0]=1.0-alB[ptr][Num_RNAType+2];
		}
	}
	
	delete [] bufA;
	delete [] bufB;
}

void Needle::GetAlignments(DIRTracks *pr, char *sA, ProbString &sB, char *alA,ProbString &alB)
{
	char *ptrA, *bufA;
	ProbChar *bufB,*ptrB;
	int i,j;
	DIRS endpt,newpt;

  i=strlen(sA);
  j=sB.Len();          
	
	int dlen=j;
	
	bufA=(char *)Malloc((i+j+1)*sizeof(char));
	bufB=new ProbChar[i+j];
	
  ptrA=bufA;
  ptrB=bufB;       
       
  while (i > 0 && j > 0)
  {
    endpt=pr[i*(dlen+1)+j].GetEntry();
    newpt=pr[i*(dlen+1)+j].GetEntryVal();
    if (endpt==DIR_DIAG)
       {
         *(ptrA++) = sA[i-1];
         *(ptrB++) = sB[j-1];
         i--;j--;                
       }
    else if (endpt==DIR_UP)
      {
        *(ptrA++) = sA[i-1];
        *(ptrB++) = gapChar;
        i--;
      }
    else if (endpt==DIR_LEFT)
      {
        *(ptrA++) = '-';
        *(ptrB++) = sB[j-1];
        j--;
      }
    else
    {	
    	cerr <<"Error No Match" << i << " " << j << endl;
    	cerr <<sA<<endl;
    	cerr <<sB.ToString()<<endl;
    	exit(0);
    }  
    pr[i*(dlen+1)+j].SetEntry(newpt);
  }
  while(i > 0)
  {
    *(ptrA++) = sA[i-1];
    *(ptrB++) = gapChar;
    i--;            
  }
  while(j > 0)
  {
    *(ptrA++) = '-';
    *(ptrB++) = sB[j - 1];
    j--;            
  }
	
	
	alB.Alloc((int) (ptrB-bufB));
	
	do {
		*(alA++)=*(--ptrA);
	} while (ptrA > bufA);		
	*alA=0;
	
	for (i=0;i<alB.Len();i++) alB[i]=*(--ptrB);

	int ptr;
	i=0;
	while (i <alB.Len() && alB[i].isGap()) alB[i++]=EndgapChar;
	j=alB.Len()-1;
	while (j>=0 && alB[j].isGap()) alB[j--]=EndgapChar;
	for (ptr=i+1;ptr<j;ptr++)
	{
		if (alB[ptr].isGap())
		{
			alB[ptr][Num_RNAType+2]=max(alB[ptr-1][Num_RNAType+2],alB[ptr][Num_RNAType+2]);
			alB[ptr][0]=1.0-alB[ptr][Num_RNAType+2];
		}
	}
	for (ptr=j-1;ptr>i;ptr--)
	{
		if (alB[ptr].isGap())
		{
			alB[ptr][Num_RNAType+2]=max(alB[ptr+1][Num_RNAType+2],alB[ptr][Num_RNAType+2]);
			alB[ptr][0]=1.0-alB[ptr][Num_RNAType+2];
		}
	}

	free(bufA);
	delete [] bufB;
}


void Needle::GetAlignments(DIRTracks *pr, ProbString &sA, ProbString &sB, ProbString &alA,ProbString &alB, int *posA, int *posB)
{
	ProbChar *bufA,*bufB,*ptrA,*ptrB;
	int i,j,pt;
	DIRS endpt,newpt;
	int count=0;
	
	i=sA.Len();
	j=sB.Len();          
	
	int dlen=j;
	
	bufA=new ProbChar[i+j];
	bufB=new ProbChar[i+j];
	
	ptrA=bufA;
	ptrB=bufB;       
       
	while (i > 0 && j > 0)
	{
		endpt=pr[i*(dlen+1)+j].GetEntry();
		newpt=pr[i*(dlen+1)+j].GetEntryVal();
		if (endpt==DIR_DIAG)
		{
			*(ptrA++) = sA[i-1];
			*(ptrB++) = sB[j-1];
			i--;j--;
			posA[i]=count;
			posB[j]=count;
		}
		else if (endpt==DIR_UP)
		{
			*(ptrA++) = sA[i-1];
			*(ptrB++) = gapChar;
			i--;
			posA[i]=count;
		}
		else if (endpt==DIR_LEFT)
		{
			*(ptrA++) = gapChar;
			*(ptrB++) = sB[j-1];
			j--;
			posB[j]=count;
		}
		else
		{	
			cerr <<"Error No Match" << i << " " << j << endl;
			cerr <<sA.ToString()<<endl;
			cerr <<sB.ToString()<<endl;
			exit(0);
		}  
		pr[i*(dlen+1)+j].SetEntry(newpt);
		count++;
	}
	while(i > 0)
	{
		*(ptrA++) = sA[i-1];
		*(ptrB++) = gapChar;
		i--;    
		posA[i]=count++;
	}
	while(j > 0)
	{
		*(ptrA++) = gapChar;
		*(ptrB++) = sB[j - 1];
		j--;            
		posB[j]=count++;
	}
	
	
	alA.Alloc((int) (ptrA-bufA));
	alB.Alloc((int) (ptrB-bufB));
	
	for (i=0;i<alA.Len();i++) 
	{
		alA[i]=*(--ptrA);
	}
	for (i=0;i<sA.Len();i++)
		posA[i]=alA.Len()-1-posA[i];
	/*	cerr << "printing out posA " <<endl;
	for (i=0;i< sA.Len();i++)
		cerr << posA[i] << " ";
	cerr <<endl;*/
	
	for (i=0;i<alB.Len();i++) 
	{
		alB[i]=*(--ptrB);
	}
	for (i=0;i<sB.Len();i++)
		posB[i]=alB.Len()-1-posB[i];
/*	cerr << "printing out posB " <<endl;
	for (i=0;i< sB.Len();i++)
		cerr << posB[i] << " ";
	cerr <<endl;*/
	
	if (ptrA!=bufA || ptrB !=bufB)
	{
		cerr << "Needle align ptr error\n"<<endl;
	}
	//convert gap to end-gap
	int ptr;
	
	i=0;
	while (i <alA.Len() && alA[i].isGap()) alA[i++]=EndgapChar;
	j=alA.Len()-1;
	while (j>=0 && alA[j].isGap()) 	alA[j--]=EndgapChar;
	for (ptr=i+1;ptr<j;ptr++)
	{
		if (alA[ptr].isGap())
		{
			alA[ptr][Num_RNAType+2]=max(alA[ptr-1][Num_RNAType+2],alA[ptr][Num_RNAType+2]);
			alA[ptr][0]=1.0-alA[ptr][Num_RNAType+2];
		}
	}
	for (ptr=j-1;ptr>i;ptr--)
	{
		if (alA[ptr].isGap())
		{
			alA[ptr][Num_RNAType+2]=max(alA[ptr+1][Num_RNAType+2],alA[ptr][Num_RNAType+2]);
			alA[ptr][0]=1.0-alA[ptr][Num_RNAType+2];
		}
	}
	
	i=0;
	while (i <alB.Len() && alB[i].isGap()) alB[i++]=EndgapChar;
	j=alB.Len()-1;
	while (j>=0 && alB[j].isGap()) alB[j--]=EndgapChar;
	for (ptr=i+1;ptr<j;ptr++)
	{
		if (alB[ptr].isGap())
		{
			alB[ptr][Num_RNAType+2]=max(alB[ptr-1][Num_RNAType+2],alB[ptr][Num_RNAType+2]);
			alB[ptr][0]=1.0-alB[ptr][Num_RNAType+2];
		}
	}
	for (ptr=j-1;ptr>i;ptr--)
	{
		if (alB[ptr].isGap())
		{
			alB[ptr][Num_RNAType+2]=max(alB[ptr+1][Num_RNAType+2],alB[ptr][Num_RNAType+2]);
			alB[ptr][0]=1.0-alB[ptr][Num_RNAType+2];
		}
	}
	
	delete [] bufA;
	delete [] bufB;
}


int NumAlign3=0;


void Needle::Align(ProbString &str1,ProbString &str2,ProbString &al1, ProbString &al2,REAL diag_fact)
{
	DIRTracks *pr;
	REAL score;
	int len1,len2;
	int i,j;
	
	NumAlign3++;
	
	len1=str1.Len();
	len2=str2.Len();
	

	pr=(DIRTracks *)Malloc((len1+1)*(len2+1)*sizeof(DIRTracks));

//	if (Diag==1)
		score=CalculateMatrix_Diag(str1, str2, len1, len2, (int) (diag_fact*max(len1,len2))+1,(int) (diag_fact*max(len1,len2))+1, pr);
//	else
//		score=CalculateMatrix(str1, str2, len1, len2, pr);

	GetAlignments(pr, str1, str2, al1, al2);

	free(pr);
}

int NumAlign1=0;
void Needle::Align(char *seq1,char *seq2,char *al1,char *al2,REAL diag_fact)
{
	DIRTracks *pr;
	REAL score;
	int len1,len2;
	int *code1, *code2;
	int i,j;

	NumAlign1++;
	
	len1=strlen(seq1);
	len2=strlen(seq2);
	code1=StrToCode(seq1);
	code2=StrToCode(seq2);
	
	pr=(DIRTracks *)Malloc((len1+1)*(len2+1)*sizeof(DIRTracks));
	

	score=CalculateMatrix_Diag(code1, code2, len1, len2, (int) (diag_fact*max(len1,len2))+1,(int) (diag_fact*max(len1,len2))+1, pr);
	GetAlignments(pr, seq1, seq2, al1, al2);

	free(pr);
	free(code1);
	free(code2);
}

int NumAlign2=0;

void Needle::Align(char *str1,ProbString &str2,char *al1, ProbString &al2,REAL diag_fact)
{
	DIRTracks *pr;
	REAL score;
	int len1,len2;
	int i,j;
	int *code1;
	
	NumAlign2++;
	
	len1=strlen(str1);
	len2=str2.Len();
	code1=StrToCode(str1);
	

	pr=(DIRTracks *)Malloc((len1+1)*(len2+1)*sizeof(DIRTracks));

	score=CalculateMatrix_Diag(code1, str2, len1, len2, (int) (diag_fact*max(len1,len2))+1,(int) (diag_fact*max(len1,len2))+1, pr);

	GetAlignments(pr, str1, str2, al1, al2);

	free(pr);
	free(code1);
}

void Needle::Align(ProbString &str1,ProbString &str2,ProbString &al1, ProbString &al2,int* pos1, int* pos2, REAL diag_fact)
{
	DIRTracks *pr;
	REAL score;
	int len1,len2;
	int i,j;
	
	NumAlign3++;
	
	len1=str1.Len();
	len2=str2.Len();
	

	pr=(DIRTracks *)Malloc((len1+1)*(len2+1)*sizeof(DIRTracks));

//	if (Diag==1)
		score=CalculateMatrix_Diag(str1, str2, len1, len2, (int) (diag_fact*max(len1,len2))+1,(int) (diag_fact*max(len1,len2))+1, pr);
//	else
//		score=CalculateMatrix(str1, str2, len1, len2, pr);

	GetAlignments(pr, str1, str2, al1, al2, pos1, pos2);

	free(pr);
}

void Needle::AlignToTemp(char *str1,ProbString &str2,char *al1,REAL diag_fact)
{
	DIRTracks *pr;
	REAL score;
	int len1,len2;
	int i,j;
	int *code1;
	ProbString al2;
	
	NumAlign2++;
	
	len1=strlen(str1);
	len2=str2.Len();
	code1=StrToCode(str1);
	

	pr=(DIRTracks *)Malloc((len1+1)*(len2+1)*sizeof(DIRTracks));

//	score=CalculateMatrix_Diag(code1, str2, len1, len2, (int) (diag_fact*max(len1,len2))+1,0, pr);
	score=CalculateMatrix_NoDgap(code1, str2, len1, len2, pr);
	
	GetAlignments(pr, str1, str2, al1, al2);

	free(pr);
	free(code1);
}
