#include <string.h>
#include <stdlib.h>
#include "util.h"
#include "ProbModel.h"
#include <iostream>
#include "global.h"

using namespace std;

REAL gaparr[]={1.0,0.0,0.0,0.0,0.0,0.0,0.0};
REAL Aarr[]={0.0,1.0,0.0,0.0,0.0,0.0,0.0};
REAL Garr[]={0.0,0.0,1.0,0.0,0.0,0.0,0.0};
REAL Carr[]={0.0,0.0,0.0,1.0,0.0,0.0,0.0};
REAL Tarr[]={0.0,0.0,0.0,0.0,1.0,0.0,0.0};
REAL Narr[]={0.0,0.0,0.0,0.0,0.0,1.0,0.0};
REAL endarr[]={0.0,0.0,0.0,0.0,0.0,0.0,1.0};

static char RNAList[]="-AGCTN_";

ProbChar gapChar(gaparr);
ProbChar EndgapChar(endarr);

ProbChar::ProbChar()
{
}

ProbChar::ProbChar(REAL arr[])
{
	Assign(arr);
}

REAL ProbChar::Dist(ProbChar &pc){
	REAL diff=probvect[0]*(1-pc[0]-pc[Num_RNAType+2]);//*0.5;
	for (int j=1; j<=Num_RNAType+1;j++)
		{
			diff+=probvect[j]*(1-pc[j]-pc[Num_RNAType+2]); //-0.5*pc[0]
		}
  return diff;
}

REAL ProbChar::Similar(ProbChar & pc){
	REAL sim=0;
	for (int j=0; j<=Num_RNAType+1;j++)
		{
			sim+=min(probvect[j],pc[j]);
		}
  return sim;
}

void ProbChar::Assign(REAL arr[])
{
	for (int i=0;i<=Num_RNAType+2;i++)
			probvect[i]=arr[i];	
}


char ProbChar::ToChar()
{
	REAL maxpr=probvect[0];
	int maxidx=0;
	for (int j=1;j<=Num_RNAType+2;j++)
		if (probvect[j]>maxpr)
		{
			maxpr=probvect[j];
			maxidx=j;
		}
	return RNAList[maxidx];
}

char ProbChar::ToWildcard()
{
	if (probvect[Num_RNAType+2] > 1.0-EPS) return '_';
	if (probvect[0] > 1.0-EPS) return '-';
	if (isGap()) return '=';
	if (probvect[Num_RNAType+2] >0.5) return ',';
	if (probvect[0]+probvect[Num_RNAType+2] >0.5) return '+';
	if (probvect[1]/(1-probvect[0]-probvect[Num_RNAType+2]) >0.5) return 'A';
	if (probvect[2]/(1-probvect[0]-probvect[Num_RNAType+2]) >0.5) return 'G';
	if (probvect[3]/(1-probvect[0]-probvect[Num_RNAType+2]) >0.5) return 'C';
	if (probvect[4]/(1-probvect[0]-probvect[Num_RNAType+2]) >0.5) return 'T';
	if (probvect[5]/(1-probvect[0]-probvect[Num_RNAType+2]) >0.5) return 'N';
	REAL max1=0;
	REAL max2=0;
	REAL max3=0;
	int idx1;
	int idx2;
	int idx3;
	int i;
	
	for (i=1;i<5;i++)
		if (probvect[i] >max1)
		{
			max1=probvect[i];
			idx1=i;
		}
	for (i=1;i<5;i++)
		if (i !=idx1 && probvect[i] >max2)
		{
			max2=probvect[i];
			idx2=i;
		}
	if ((max1+max2)/(1-probvect[0]-probvect[Num_RNAType+2]) >0.5)
	{
		if (idx1==1 && idx2==2 || idx1==2 && idx2==1) return 'R';
		if (idx1==1 && idx2==3 || idx1==3 && idx2==1) return 'M';
		if (idx1==1 && idx2==4 || idx1==4 && idx2==1) return 'W';
		if (idx1==2 && idx2==3 || idx1==3 && idx2==2) return 'S';
		if (idx1==2 && idx2==4 || idx1==4 && idx2==2) return 'K';
		if (idx1==3 && idx2==4 || idx1==4 && idx2==3) return 'Y';
	}
	
	for (i=1;i<5;i++)
		if (i !=idx1 && i!=idx2 && probvect[i] >max3)
		{
			max3=probvect[i];
			idx3=i;
		}
	if ((max1+max2+max3)/(1-probvect[0]-probvect[Num_RNAType+2]) >0.75)
	{
		if (10-idx1-idx2-idx3 ==1) return 'B';
		if (10-idx1-idx2-idx3 ==2) return 'D';
		if (10-idx1-idx2-idx3 ==3) return 'H';
		if (10-idx1-idx2-idx3 ==4) return 'V';
	}
	return 'N';
}

void ProbString::Clear()
{
		if (seq !=NULL) 
		{			
			delete [] seq;
			seq=NULL;
		}
}

void ProbString::Alloc(int len)
{
	Clear();
	this->len=len;
	seq=new ProbChar[len];
}

ProbString::~ProbString()
{
	Clear();
}

ProbString & ProbString::operator =(ProbString &ps)
{
	Alloc(ps.Len());
	
	for (int i=0;i<len;i++)
		seq[i]=ps[i];
	return ps;
}

ProbString::ProbString(char *strseq)
{
	seq=NULL;
	Alloc(strlen(strseq));
	for (int i=0;i<len;i++)
	{
		switch (strseq[i]){
			case '-': seq[i].Assign(gaparr);
				break;
			case 'A': seq[i].Assign(Aarr);
				break;
			case 'G': seq[i].Assign(Garr);
				break;
			case 'C': seq[i].Assign(Carr);
				break;
			case 'T': seq[i].Assign(Tarr);
				break;
			case 'N': seq[i].Assign(Narr);
				break;			
			default:
				break;				
		}
	}
		int ptr=0;
	while (ptr <len && seq[ptr].isGap())
	{
		seq[ptr++].Assign(endarr);
	}
	ptr=len-1;
	while (ptr >=0 && seq[ptr].isGap())
	{
		seq[ptr--].Assign(endarr);
	}
}

void ProbString::FromString(char *strseq)
{
	Alloc(strlen(strseq));
	for (int i=0;i<len;i++)
	{
		switch (strseq[i]){
			case '-': seq[i].Assign(gaparr);
				break;
			case 'A': seq[i].Assign(Aarr);
				break;
			case 'G': seq[i].Assign(Garr);
				break;
			case 'C': seq[i].Assign(Carr);
				break;
			case 'T': seq[i].Assign(Tarr);
				break;
			case 'N': seq[i].Assign(Narr);
				break;			
			default:
				break;				
		}
	}
	int ptr=0;
	while (ptr <len && seq[ptr].isGap())
	{
		seq[ptr++].Assign(endarr);
	}
	ptr=len-1;
	while (ptr >=0 && seq[ptr].isGap())
	{
		seq[ptr--].Assign(endarr);
	}
}

REAL ProbString::AveDist(ProbString &str)
{
	int gap_alert1 = 0;
	int gap_alert2 = 0;

	int start_seq1 = 0;
	int start_seq2 = 0;
	int end_seq1 = this->len;
	int end_seq2 = str.Len();
	int k;
	REAL residuecount, distance;
	REAL accgap=0.0;
	REAL prob_res;
	
	residuecount = distance = 0.0;

	k = 0;
	while(seq[k].isGap()) k++;
	start_seq1 =  k;
	k = 0;
	while(str[k].isGap()) k++;
	start_seq2 = k;
	
	k = end_seq1-1;
	while(seq[k].isGap()) k--;
	end_seq1 = k+1; 
	  	
	k = end_seq2-1;
	while(str[k].isGap()) k--;
	end_seq2 = k+1; 

	gap_alert1 = 0;
	gap_alert2 = 0;
	
	for( k=max(start_seq1, start_seq2); k <min(end_seq1, end_seq2); k++) 
	{

		if (seq[k].isGap() && str[k].isGap())
		{
			continue;
		}

		if(seq[k].isUnknown() && str[k].isUnknown())
		{
			if (gap_alert1>0) 
			{
				residuecount+=prob_res;
				//residuecount+=1;
				distance+=accgap/gap_alert1;
				accgap=0.0;
				prob_res=1.0;
			}
			if (gap_alert2>0) 
			{
				residuecount+=prob_res;
				//residuecount+=1;
				distance+=accgap/gap_alert2;
				accgap=0.0;
				prob_res=1.0;
			}
			residuecount+=1.0;
			gap_alert1 = 0;
			gap_alert2 = 0;
			continue;
		}

		if (seq[k].isGap())
		{
			if(gap_alert1 == 0)
			{
				if (gap_alert2>0) 
				{
					residuecount+=prob_res;
					//residuecount+=1;
					distance+=accgap/gap_alert2;
				}
				accgap=0.0;
				gap_alert2 = 0;
				prob_res=0.0;
			}
			accgap+=seq[k].Dist(str[k]);
			prob_res=max(prob_res,(1-str[k].ProbEnd())*(1-seq[k].ProbEnd()));
			gap_alert1++;
			continue;
		}

		if ( str[k].isGap())
		{
			if(gap_alert2 == 0)
			{
				if (gap_alert1>0) 
				{
					residuecount+=prob_res;
					//residuecount+=1;
					distance+=accgap/gap_alert1;
				}
				prob_res=0.0;
				accgap=0.0;
				gap_alert1 = 0;
			}
			accgap+=seq[k].Dist(str[k]);
			prob_res=max(prob_res,(1-str[k].ProbEnd())*(1-seq[k].ProbEnd()));
			gap_alert2++;
			continue;
		}

		distance += seq[k].Dist(str[k]);
		residuecount+=(1.0-seq[k].ProbEnd())*(1.0-str[k].ProbEnd());
		if (gap_alert1>0) 
		{
			distance+=accgap/gap_alert1;
			residuecount+=prob_res;
			//residuecount+=1;
			accgap=0.0;
			prob_res=1.0;
		}
		if (gap_alert2>0) 
		{
			distance+=accgap/gap_alert2;
			residuecount+=prob_res;
			//residuecount+=1;
			accgap=0.0;
			prob_res=1.0;
		}
		gap_alert1 = 0;
		gap_alert2 = 0;
	}

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

REAL ProbString::AveDist(char *str)
{
	int gap_alert1 = 0;
	int gap_alert2 = 0;

	int start_seq1 = 0;
	int start_seq2 = 0;
	int end_seq1 = this->len;
	int end_seq2 = strlen(str);
	int k;
	REAL residuecount, distance;
	REAL accgap=0.0;
	REAL prob_res;
	
	residuecount = distance = 0.0;

	k = 0;
	while(seq[k].isGap()) k++;
	start_seq1 =  k;
	k = 0;
	while(str[k]=='-') k++;
	start_seq2 = k;
	k = end_seq1-1;
	while(seq[k].isGap()) k--;
	end_seq1 = k+1; 
	k = end_seq2-1;
	while(str[k]=='-') k--;
	end_seq2 = k+1; 

	gap_alert1 = 0;
	gap_alert2 = 0;

	for( k=max(start_seq1, start_seq2); k <min(end_seq1, end_seq2); k++) {

		if (seq[k].isGap() && str[k]=='-')
		{
			continue;
		}

	  	if(seq[k].isUnknown() && str[k]=='N')
		{
			if (gap_alert1>0) 
			{	
				distance+=accgap/gap_alert1;
				residuecount+=prob_res;
				//residuecount+=1;
				accgap=0.0;
				prob_res=1.0;
			}
			if (gap_alert2>0) 
			{
				distance+=accgap/gap_alert2;
				//residuecount+=prob_res;
				residuecount+=1;
				accgap=0.0;
				prob_res=1.0;
			}
			residuecount+=1;
			gap_alert1 = 0;
			gap_alert2 = 0;
			continue;
		}

		if ( seq[k].isGap())
		{
			if(gap_alert1 == 0)
			{
				if (gap_alert2>0) 
				{	
					distance+=accgap/gap_alert2;
					residuecount+=prob_res;
					//residuecount+=1;
				}
				accgap=0.0;
				prob_res=0;
				gap_alert2 = 0;
			}
			accgap+=1-seq[k].ProbEnd();
			prob_res=max(prob_res,1-seq[k].ProbEnd());
			gap_alert1++;
			continue;
		}

		
    	if ( str[k]=='-')
		{
			if(gap_alert2 == 0)
			{
				if (gap_alert1 >0)
				{
					distance+=accgap/gap_alert1;
					residuecount+=prob_res;
					//residuecount+=1;
				}
				accgap=0.0;
				prob_res=0.0;
				gap_alert1 = 0;
			}
           	accgap+=1-seq[k].ProbGap();
			prob_res=max(prob_res,1-seq[k].ProbEnd());
			gap_alert2++;
			continue;
		}

   	 	int idx=(int) (strchr(RNAList,str[k])-RNAList);
     	distance+=1-seq[k][idx]-seq[k].ProbEnd();
		residuecount+=1.0-seq[k].ProbEnd();
		if (gap_alert1>0) 
		{
			distance+=accgap/gap_alert1;
			residuecount+=prob_res;
			//residuecount+=1;
			accgap=0.0;
			prob_res=1.0;
		}
		if (gap_alert2>0) 
		{	
			distance+=accgap/gap_alert2;
			residuecount+=prob_res;
			//residuecount+=1;
			accgap=0.0;
			prob_res=1.0;
		}
		gap_alert1 = 0;
		gap_alert2 = 0;
	}

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



int CompRealD(const void * a, const void * b)
{
	if (*(REAL *)a > *(REAL *)b) return -1;
	if (*(REAL *)a < *(REAL *)b) return 1;
	return 0;
}


int ProbString::findcutofflen(REAL *difvec,int difcnt,REAL startprob,REAL pval)
{
	REAL prob=startprob;
	int difptr=0; //max(difstart,0);
	REAL comb;
	int i,j,k;
	int ret;
	bool moveon;
		
	while (prob >pval & difptr<difcnt)
	{
		prob *=difvec[difptr++];
	}	
	difptr--;
	prob/=difvec[difptr];
		
	if (difptr < difcnt-1)
	{
		comb=0.0;
		for (i=difptr;i<difcnt;i++)
			comb+=difvec[i];
		while (comb*prob >pval & difptr <difcnt)
		{
			prob *=difvec[difptr++];
			comb=0.0;
			for (i=difptr;i<difcnt;i++)
				comb+=difvec[i];
		}
	}	
	ret=difptr;
	if (difptr < difcnt-2)
	{
		comb=0.0;
		for (i=difptr;i<difcnt;i++)
		{
			for (j=i+1;j<difcnt;j++)
				comb+=difvec[i]*difvec[j];
		}	
		moveon=false;
		while (comb*prob >pval & difptr <difcnt)
		{
			moveon=true;
			prob *=difvec[difptr++];
			comb=0.0;
			for (i=difptr;i<difcnt;i++)
			{
				for (j=i+1;j<difcnt;j++)
					comb+=difvec[i]*difvec[j];
			}	
		}
	}	
	if (moveon) ret=difptr+1;
	
	/*if (difptr < difcnt-3)
	{
		comb=0.0;
		for (i=difptr;i<difcnt;i++)
		{
			for (j=i+1;j<difcnt;j++)
				for (k=j+1;k<difcnt;k++)
					comb+=difvec[i]*difvec[j]*difvec[k];
		}	
		moveon=false;
		while (comb*prob >pval & difptr <difcnt)
		{
			moveon=true;
			prob *=difvec[difptr++];
			comb=0.0;
			for (i=difptr;i<difcnt;i++)
			{
				for (j=i+1;j<difcnt;j++)
					for (k=j+1;k<difcnt;k++)
						comb+=difvec[i]*difvec[j]*difvec[k];
			}	
		}
	}	
	if (moveon) ret=difptr+2;*/
	
//	return max(ret,difstart);			
	return ret;
}

REAL ProbString::ProbMaxDist(ProbString &str, REAL pval)
{
	int gap_alert1 = 0;
	int gap_alert2 = 0;

	int start_seq1 = 0;
	int start_seq2 = 0;
	int end_seq1 = this->len;
	int end_seq2 = str.Len();
	int k;
	REAL residuecount, distance;
	REAL accgap=0.0;
	REAL startprob=1.0;

	REAL *difvec=(REAL *)Malloc(max(end_seq1,end_seq2)*sizeof(REAL));
	int difcnt=0;
	residuecount = distance = 0.0;

	k = 0;
	while(seq[k].isGap()) k++;
	start_seq1 =  k;
	k = 0;
	while(str[k].isGap()) k++;
	start_seq2 = k;
	k = end_seq1-1;
	while(seq[k].isGap()) k--;
	end_seq1 = k+1; 
	k = end_seq2-1;
	while(str[k].isGap()) k--;
	end_seq2 = k+1; 
	gap_alert1 = 0;
	gap_alert2 = 0;

	for( k=max(start_seq1, start_seq2); k <min(end_seq1, end_seq2); k++) {

		if (seq[k].isGap() && str[k].isGap())
		{
			continue;
		}

	  	if(seq[k].isUnknown() && str[k].isUnknown())
		{
			if (gap_alert1>0) 
			{
				distance+=accgap/gap_alert1;
				startprob=min(startprob,accgap/gap_alert1);
				accgap=0.0;
			}
			if (gap_alert2>0) 
			{
				distance+=accgap/gap_alert2;
				startprob=min(startprob,accgap/gap_alert2);
				accgap=0.0;
			}
			gap_alert1 = 0;
			gap_alert2 = 0;
			continue;
		}

    	if ( seq[k].isGap())
		{
			if(gap_alert1 == 0)
			{
				residuecount+=1.0;
				if (gap_alert2>0) 
				{
					distance+=accgap/gap_alert2;
					startprob=min(startprob,accgap/gap_alert2);
					accgap=0.0;
				}
				gap_alert2 = 0;
			}
			accgap+=1-str[k][0];
					//startprob=min(startprob,1-str[k][0]);
					//difvec[difcnt++]=1-str[k][0];
			gap_alert1++;
			continue;
		}

    	if ( str[k].isGap())
		{
			if(gap_alert2 == 0)
			{
				residuecount+=1.0;
				if (gap_alert1>0) 
				{
					distance+=accgap/gap_alert1;
					startprob=min(startprob,accgap/gap_alert1);
					accgap=0.0;
				}
				gap_alert1 = 0;
			}
			gap_alert2 = 1;
			accgap+=1-seq[k][0];
					//startprob=min(startprob,1-seq[k][0]);
					//difvec[difcnt++]=1-seq[k][0];
			continue;
		}
		
    	 //difvec[difcnt++]= 1-seq[k].Similar(str[k]);
    	difvec[difcnt++]= seq[k].Dist(str[k]);
        residuecount+=1.0;
		if (gap_alert1>0) 
		{
			distance+=accgap/gap_alert1;
			startprob=min(startprob,accgap/gap_alert1);
			accgap=0.0;
		}
		if (gap_alert2>0) 
		{
			distance+=accgap/gap_alert2;
			startprob=min(startprob,accgap/gap_alert2);
			accgap=0.0;
		}

		gap_alert1 = 0;
		gap_alert2 = 0;
	}

  if (residuecount > 0) {
		qsort(difvec,difcnt,sizeof(REAL),CompRealD);

		distance +=findcutofflen(difvec,difcnt,startprob,pval);
	
		distance = distance / residuecount;
  }
  else {
		distance = 1.0;
  }
 	free(difvec);
  return distance;  
}

REAL ProbString::ProbMaxDist(char *str, REAL pval)
{

	int gap_alert1 = 0;
	int gap_alert2 = 0;

	int start_seq1 = 0;
	int start_seq2 = 0;
	int end_seq1 = this->len;
	int end_seq2 = strlen(str);
	int k;
	REAL residuecount, distance;
	REAL startprob=1.0;
	REAL accgap=0.0;
	REAL *difvec=(REAL *)Malloc(max(end_seq1,end_seq2)*sizeof(REAL));
	int difcnt=0;

	residuecount = distance = 0.0;

	k = 0;
	while(seq[k].isGap()) k++;
	start_seq1 =  k;
	
	k = 0;
	while(str[k]=='-') k++;
	start_seq2 = k;
	
	k = end_seq1-1;
	while(seq[k].isGap()) k--;
	end_seq1 = k+1; 
	  	
	k = end_seq2-1;
	while(str[k]=='-') k--;
	end_seq2 = k+1; 
	
	gap_alert1 = 0;
	gap_alert2 = 0;

	for( k=max(start_seq1, start_seq2); k <min(end_seq1, end_seq2); k++) {

		if (seq[k].isGap() && str[k]=='-')
		{
			continue;
		}

	  	if(seq[k].isUnknown() && str[k]=='N')
		{
			if (gap_alert2>0) 
			{
				distance+=accgap/gap_alert2;
				startprob=min(startprob,accgap/gap_alert2);
				accgap=0.0;
			}
			gap_alert1 = 0;
			gap_alert2 = 0;
			continue;
		}

    	if ( seq[k].isGap())
		{
			if(gap_alert1 == 0)
			{
				residuecount+=1.0;
				if (gap_alert2>0) 
				{
					distance+=accgap/gap_alert2;
					startprob=min(startprob,accgap/gap_alert2);
					accgap=0.0;
				}
				distance+=1;
				gap_alert2 = 0;
				//difvec[difcnt++]=1.0;
			}
			gap_alert1 = 1;
			continue;
		}

    	if ( str[k]=='-')
		{
			if(gap_alert2 == 0)
			{
				residuecount+=1.0;
				gap_alert1 = 0;
			}
			accgap+=1-seq[k][0];
					//startprob=min(startprob,1-seq[k][0]);
					//difvec[difcnt++]=1-seq[k][0];
			gap_alert2++;
			continue;
		}
		
    	 //difvec[difcnt++]= 1-seq[k].Similar(str[k]);
   	 	int idx=(int) (strchr(RNAList,str[k])-RNAList);
    	difvec[difcnt++]= 1-seq[k][idx];
		residuecount+=1.0;
		if (gap_alert2>0) 
		{
			distance+=accgap/gap_alert2;
			startprob=min(startprob,accgap/gap_alert2);
			accgap=0.0;
		}
		gap_alert1 = 0;
		gap_alert2 = 0;
	}

  if (residuecount > 0) {
		qsort(difvec,difcnt,sizeof(REAL),CompRealD);

		distance +=findcutofflen(difvec,difcnt,startprob,pval);
		distance = distance / residuecount;
  }
  else {
		distance = 1.0;
  }
 	free(difvec);
  return distance;  
}

void ProbString::operator +=(ProbString &ps)
{
	if (len !=ps.Len())
	{
		cerr << "Cannot combine ProbString of different length" <<endl;
		return;
	}
	for (int i=0;i<len;i++)
		seq[i]+=ps[i];
} 
	
void ProbString::operator *=(REAL w)
{
	for (int i=0;i<len;i++)
		seq[i]*=w;	
} 

void ProbString::Add(char *str, REAL wt)
{
	int idx;
	for (int i=0;i<len;i++)
	{
		idx=(int) (strchr(RNAList,str[i])-RNAList);
		seq[i][idx]+=wt;
	}
}

REAL ProbString::MinDist(ProbString &str){
	return 0;
}

bool ProbString::Rectify(REAL thres)
{
	int numreserve=0;
	int i,j,ptr;
	for (i=0;i<len;i++)
	{
		int cntlive=0;
		for (j=1; j<=Num_RNAType+1;j++)
		{
			if (seq[i][j] > thres) cntlive++;
		}
		if (cntlive>0) numreserve++;
	}
	ProbChar *newseq = new ProbChar[numreserve];

	ptr=0;
	for (i=0;i<len;i++)
	{
		REAL denorm=0;
		int cntlive=0;
		for (j=1; j<=Num_RNAType+1;j++)
		{
			if (seq[i][j] > thres) cntlive++;
		}
		
		for (j=0; j<=Num_RNAType+2;j++)
		{
			if (seq[i][j] > thres) // || (j>0 && cntlive>1 &&seq[i][j] > thres/2))
			{
				newseq[ptr][j]=seq[i][j];
				denorm+=seq[i][j];
			}
			else
			{
				newseq[ptr][j]=0.0;
			}
		}
		if (cntlive>0)
		{
			newseq[ptr] *=1.0/denorm;
			ptr++;
			if (ptr >=numreserve) break;
		}
	}
	Clear();
	seq=newseq;
	bool dorec=(numreserve <len);
	len=numreserve;
	return dorec;
}


int ProbString::DOF()
{
	int res=0;
	for (int i=0;i<len;i++)
	{
		int cnt=0;
		for (int j=0;j<=Num_RNAType;j++)
		{
			if (seq[i][j]>MIN_PROB) cnt++;
		}
		if (cnt>1) res++;
	}
	return res;
}

char* ProbString::ToString()
{
	char *buff;
	buff=(char *)Malloc((len+1)*sizeof(char));
	for (int i=0;i<len;i++)
	{
		buff[i]=seq[i].ToChar();
	}
	buff[len]='\0';
	return buff;
}

char* ProbString::ToWildcard()
{
	char *buff;
	buff=(char *)Malloc((len+1)*sizeof(char));
	for (int i=0;i<len;i++)
	{
		buff[i]=seq[i].ToWildcard();
	}
	buff[len]='\0';
	return buff;
}

fstream &operator <<(fstream &out, ProbChar &chr)
{
	out << (int) ((chr[0]+chr[Num_RNAType+2]+1/510)*255) << " ";
	
	for (int i=1;i<=Num_RNAType;i++)
	{
		out << (int) ((chr[i]+1/510)*255) << " ";
	}
	return out;
}

fstream &operator <<(fstream &out, ProbString &str)
{
	out << str.Len() << " ";
	for (int i=0;i<str.Len(); i++)
 			out<<str[i];
 	out << endl;
 return out;
}

fstream &operator >>(fstream &in,ProbChar &chr)
{
	int ntprob;
	for (int i=0;i<=Num_RNAType;i++)
	{
		in >>  ntprob;
		chr[i]=(REAL) ntprob/255;
	}
	chr[Num_RNAType+1]=0;
	chr[Num_RNAType+2]=0;
	return in;
	
}

fstream &operator >>(fstream &in,ProbString &str)
{
	int len,i;
	in >> len;
	str.Alloc(len);
	for (i=0;i<len;i++)
	{
		in >>str[i];
	}
	return in;
}

