#include "kmer.h"
#include "ProbCluster.h"
#include "needle.h"
#include <iostream>
#include <fstream>
#include <time.h>
#include <math.h>

#include "global.h"
#include "MinHeap.h"
using namespace std;

extern bool logon;

extern int *SeqID;

extern int NumAlign1;
extern int NumAlign2;
extern int NumAlign3;
extern long int Kcomp;
extern long int Kcomp2;

extern double AveSeqLen;

ProbCluster::ProbCluster(REAL lmin,REAL lmax)
{
	root=NULL;
	levelinc=Global::level_inc;
	levelmax=lmax*levelinc;
	levelmin=levelmax;
	
	if (Global::UseApprox)
		while (levelmin > min(lmin,(REAL)(3.0/AveSeqLen))) levelmin=levelmin/levelinc;
	else
		while (levelmin > min(lmin/levelinc,(REAL)(3.0/AveSeqLen))) levelmin=levelmin/levelinc;
	
	maxdist= Global::level_max*Global::level_inc;
}

clock_t lastt,thist;

void ProbCluster::InitCluster(int SeqNum)
{
	clRec=new ClusterRecord[2*SeqNum];
	clistTop=SeqNum;
	NumCls=SeqNum;

}

void ProbCluster::Finish()
{
	delete [] clRec;
}

void ProbCluster::AddRecord(unsigned int uid, Tree *node)
{
	clRec[uid].Node=node;
	clRec[uid].NumSeqs=node->NumSeqs();
	clRec[uid].NNdist=maxdist;
	clRec[uid].clsid=-1;
	clRec[uid].NNseq=-1;
	clRec[uid].ID=uid;
	clRec[uid].distErr=0.0;
	clRec[uid].tagscan=0;
}

Tree * ProbCluster::CreateBranch(unsigned int uid,int frq,REAL levelup)
{
	Tree *base=new Tree(uid,frq);
	Tree *upper;
	REAL level=levelmin;
	while (level <=levelup*1.001)
	{
		upper=new Tree(uid,frq);
		upper->AddChild(base,false);
		//upper->NumSeqs()=frq;
		upper->SingleSeq=1;
		
		upper->SetThres(level);
		level*=levelinc;
		base=upper;
	}
	return base;
}	

Tree * ProbCluster::CreateBranch(Tree *node,REAL levelup)
{
	Tree *base=node;
	Tree *upper;
	REAL level=levelmin;
	while (level <=levelup*1.001)
	{
		upper=new Tree;
		upper->AddChild(base);
		if (node->SingleSeq)
		{
			upper->SingleSeq=1;
		}
		else
		{
			upper->Center()=base->Center();
			upper->SingleSeq=0;
		}
		upper->UID=node->UID;

		upper->SetThres(level);
		level*=levelinc;
		base=upper;
	}
		
	return base;
}	

void ProbCluster::AddSeq(unsigned int uid,int frq)
{
	
	Tree *current,*bestchild;
	REAL nrdist;
	ProbString al;
	if (root==NULL)
	{
		root=CreateBranch(uid,frq,levelmax*levelinc);	
		current=root;
		while (!current->IsLeaf()) 
		{
			current=current->FirstChild();
		}
		AddRecord(uid,current);
		clRec[uid].Seqlist.clear();
		clRec[uid].Seqlist.insert(uid);
		return;
	}
	
	current=root;
//	current->NumSeqs()+=frq;

	do{
		if (current->IsLeaf()) cerr <<"Shouldn't Get to Leaf" <<endl;
		
		if (current->BottomLevel())
		{
			bestchild=new Tree(uid,frq);
			current->AddChild(bestchild);
			break;
		}
		
		if (current->NumChildren()==1)
		{
			bestchild=current->FirstChild();
			if (current==root)
			{
				if (current->SingleSeq)
				{
					Global::needle->Align(SeqStrs[uid],SeqStrs[bestchild->UID],Global::buf1,Global::buf2,Global::DiagRate);
					nrdist=calcDistance(Global::buf1,Global::buf2);	
				}
				else
				{
					Global::needle->Align(SeqStrs[uid],bestchild->Center(),Global::buf1,al,Global::DiagRate);
					if (Global::Dist_Style==DS_Ave)	
						nrdist=al.AveDist(Global::buf1);
					else
						nrdist=al.ProbMaxDist(Global::buf1,0.1/bestchild->NumSeqs());
				}
				Global::cache->Store(uid,bestchild->UID,nrdist);
			}
		}
		else
		{
			nrdist=current->FindSpanChild(uid,bestchild);
		}

		
		if (bestchild !=NULL && nrdist < bestchild->GetThres())
		{
			current->NumSeqs()+=frq;
			current=bestchild;
		}
		else
		{
			bestchild=CreateBranch(uid,frq,current->GetThres()/levelinc);
			current->AddChild(bestchild);
			break;
		}
	}while(1);
		
	while (!bestchild->IsLeaf()) 
	{
		bestchild=bestchild->FirstChild();
	}
	AddRecord(uid,bestchild);
	clRec[uid].Seqlist.clear();
	clRec[uid].Seqlist.insert(uid);
	return;
}


void ProbCluster::AddNode(Tree *node,Tree *top,int id1,int id2,int sq1,int sq2)
{
	Tree *bestchild;
	REAL nrdist;
	
	ProbString al1,al2;

	if (top->NumChildren()==0)
	{
		bestchild=CreateBranch(node,top->GetThres()/levelinc);
		top->AddChild(bestchild,false);
		top->UID=node->UID;
		if (!node->SingleSeq)
			{
				top->SingleSeq=0;
				top->Center()=node->Center();
			}
			else
				top->SingleSeq=1;
		return;
	}

	Tree *current=top;

	do{
		if (current->IsLeaf()) cerr <<"Shouldn't Get to Leaf" <<endl;
		if (current->BottomLevel())
		{
			current->AddChild(node,(current !=top));
			break;
		}
		
		if (current->NumChildren()==1)
		{
			bestchild=current->FirstChild();
			if (current==top )
			{
				if (!QuickDist(id1,id2,current->UID,node->UID,sq1,sq2,nrdist))
					nrdist=current->NeedleDist(node);
			}
		}
		else
		{
			if (node->SingleSeq)
				nrdist=current->FindSpanChild(node->UID,bestchild);
			else
				nrdist=current->FindSpanChild(node,bestchild,id1,id2,sq1,sq2);
		}
		
		if (bestchild !=NULL && nrdist < bestchild->GetThres())
		{
			if (current->NumChildren()==1)
			{
				if (current !=top)
					current->NumSeqs()+=node->NumSeqs();
				bestchild->NumSeqs()=current->NumSeqs();
				current=bestchild;
			}
			else
			{
				current=bestchild;
				if (current !=top)
					current->NumSeqs()+=node->NumSeqs();
			}
		}
		else
		{
			bestchild=CreateBranch(node,current->GetThres()/levelinc);
			current->AddChild(bestchild,(current !=top));
			break;
		}
	}while(1);
}


Tree *ProbCluster::MergeNode(Tree *node1,Tree *node2, unsigned int uid)
{
	ProbString pr1,pr2;
	
	if (node1->SingleSeq) 
		pr1.FromString(SeqStrs[node1->UID]); 
	else 
		pr1=node1->Center();
	if (node2->SingleSeq) 
		pr2.FromString(SeqStrs[node2->UID]); 
	else 
		pr2=node2->Center();

	ProbString al1,al2;
	if (Global::ShowAlign)
	{
		int* pos1,*pos2;
		int i;
		set <unsigned int>::iterator it;
		
		pos1=(int *)Malloc(pr1.Len()*sizeof(int));
		pos2=(int *)Malloc(pr2.Len()*sizeof(int));

		Global::needle->Align(pr1,pr2,al1,al2,pos1,pos2,Global::DiagRate);
		
		char *temp=(char *)Malloc((al1.Len()+1)*sizeof(char));
		
		for (it=clRec[node1->UID].Seqlist.begin();it!=clRec[node1->UID].Seqlist.end();it++)
		{
			if (strlen(AlignedSeqs[SeqID[*it]]) != pr1.Len())
			{
				cerr << "Len not match 1 " << strlen(AlignedSeqs[SeqID[*it]]) << " " << pr1.Len() << endl;
			}
			for (i=0;i<al1.Len();i++) temp[i]='-';
			temp[al1.Len()]=0;
			for (i=0;i<pr1.Len();i++) 
			{
				temp[pos1[i]]=AlignedSeqs[SeqID[*it]][i];
			}
			free(AlignedSeqs[SeqID[*it]]);
			AlignedSeqs[SeqID[*it]]=(char *)Malloc((al1.Len()+1)*sizeof(char));
			strcpy(AlignedSeqs[SeqID[*it]],temp);
		}	

		for (it=clRec[node2->UID].Seqlist.begin();it!=clRec[node2->UID].Seqlist.end();it++)
		{
			if (strlen(AlignedSeqs[SeqID[*it]]) != pr2.Len())
			{
				cerr << "Len not match 2 " << strlen(AlignedSeqs[SeqID[*it]]) << " " << pr2.Len() <<endl;
			}
			for (i=0;i<al1.Len();i++) temp[i]='-';
			temp[al1.Len()]=0;
			for (i=0;i<pr2.Len();i++) 
			{
				temp[pos2[i]]=AlignedSeqs[SeqID[*it]][i];
			}
			free(AlignedSeqs[SeqID[*it]]);
			AlignedSeqs[SeqID[*it]]=(char *)Malloc((al1.Len()+1)*sizeof(char));
			strcpy(AlignedSeqs[SeqID[*it]],temp);
		}	
		
		free(pos1);
		free(pos2);
	}
	else
	{
		Global::needle->Align(pr1,pr2,al1,al2,Global::DiagRate);
	}
	
	Tree *res=new Tree();
	res->NumSeqs()=node1->NumSeqs();
	res->AppendProbVect(al1,al2,node2->NumSeqs());
	res->NumSeqs()+=node2->NumSeqs();
	res->UID=uid;
	//res->BuildKmerCode();	
	bool dorec=false;
	if (res->NumSeqs() > 1.0/Global::ErrorRate) 
	{
		dorec=res->Center().Rectify(Global::ErrorRate);
	}
	
/*	if (dorec)
	{
		res->BuildKmerCode();
	}
	else*/
	{
		if (node1->SingleSeq)
		{
			if (node2->SingleSeq)
				KmerCodes[uid-SeqNum]=Global::kmer->KmerAdd(KmerTabs[node1->UID],KmerTabs[node2->UID],node1->NumSeqs(),node2->NumSeqs());
			else
				KmerCodes[uid-SeqNum]=Global::kmer->KmerAdd(KmerTabs[node1->UID],KmerCodes[node2->UID - SeqNum],node1->NumSeqs(),node2->NumSeqs());
		}
		else
		{
			if (node2->SingleSeq)
				KmerCodes[uid-SeqNum]=Global::kmer->KmerAdd(KmerTabs[node2->UID],KmerCodes[node1->UID-SeqNum],node2->NumSeqs(),node1->NumSeqs());
			else
				KmerCodes[uid-SeqNum]=Global::kmer->KmerAdd(KmerCodes[node1->UID - SeqNum],KmerCodes[node2->UID - SeqNum],node1->NumSeqs(),node2->NumSeqs());
		}
	}
	res->SingleSeq=0;
	return res;
}

Tree *ProbCluster::FindCommonParent(Tree *node1,Tree *node2)
{
	Tree *p1=node1->GetParent();
	Tree *p2=node2->GetParent();
	if (p1==p2) return p1;
	while (p1!=root && p1->GetThres() < p2->GetThres()*0.999) p1=p1->GetParent(); 
	while (p2!=root && p2->GetThres() < p1->GetThres()*0.999) p2=p2->GetParent(); 
	while (p1!=root && p2!=root && p1!=p2)
	{
		p1=p1->GetParent();
		p2=p2->GetParent();
	}
	if (p1!=p2) return root; else return p1;
}

void ProbCluster::DeleteNode(Tree *node,Tree *top)
{
	Tree *par=node->GetParent();
	int frq=node->NumSeqs();
	unsigned int ui=node->UID;
	
	while (par !=top && par->NumChildren()==1)
	{
		par->DeleteChild(node);
		node=par;
		par=par->GetParent();
	}
	par->DeleteChild(node);
	
	if (par->UID !=ui)
	{
		if (ui <SeqNum)
		{
			free(KmerTabs[ui]);
			free(KmerSeqs[ui]);
			KmerTabs[ui]=NULL;
			KmerSeqs[ui]=NULL;
		}
		else
		{
			free(KmerCodes[ui - SeqNum]);
			KmerCodes[ui - SeqNum]=NULL;
		}
	}
	
	while (par !=top)
	{
		par->NumSeqs()-=frq;
		par=par->GetParent();
	}
	
}


void ProbCluster::CondenseBottom()
{
	vector<Tree *> tvec;
	vector<Tree *>::iterator it;
	vector<Tree *> tchild;
	vector<Tree *>::iterator it2;

	root->ListChildrenAt(levelmin,tvec);
	
	REAL aveFreq=(SumFreq+0.0)/tvec.size();
	int NumCon=0;
	for (it=tvec.begin();it<tvec.end();it++)
	{
		Tree *cur=*it;
		//if (cur->NumSeqs() > 3*aveFreq && cur->NumChildren() >10)
		if (cur->NumChildren() > 1)
		{
			NumCon++;
			AddRecord(clistTop,cur);
			tchild.clear();
			cur->ListLeaf(tchild);
			int NumLeaf=tchild.size();
			
			clRec[clistTop].Seqlist.clear();
			for (it2=tchild.begin();it2!=tchild.end();it2++)
			{
				Tree *leaf=*it2;
				clRec[clistTop].Seqlist.insert(leaf->UID);
				clRec[leaf->UID].Seqlist.clear();
				clRec[leaf->UID].clsid=clistTop;
			}
			cur->Collapse(clistTop);
			clRec[clistTop].Node=cur->FirstChild();
			clistTop++;
			NumCls-=(NumLeaf-1);
			
			if (Global::ShowAlign)
			{
				char *temp;
				set <unsigned int>::iterator it3;
				for (it3=clRec[clistTop-1].Seqlist.begin();it3!=clRec[clistTop-1].Seqlist.end();it3++)
				{
					temp=(char *)Malloc(1000*sizeof(char));
					Global::needle->AlignToTemp(SeqStrs[*it3],cur->FirstChild()->Center(),temp,Global::DiagRate);
					for (int i=strlen(temp); i<cur->Center().Len(); i++)
						temp[i]='-';
					temp[cur->FirstChild()->Center().Len()]=0;
					
					free(AlignedSeqs[SeqID[*it3]]);
					AlignedSeqs[SeqID[*it3]]=temp;
				}
			}
		}
	}
}


void ProbCluster::NNSearch(int uid, bool exactsearch, int id1,int id2,int sq1,int sq2)
{

	if (clRec[uid].clsid >=0) return;
	vector<Tree *> tvec;
	Tree *thisnode=clRec[uid].Node;
	MinHeap heap;
	REAL dist;
	Tree *top;
	REAL mindist=min(maxdist,clRec[uid].NNdist);
	int NNID;
	if (clRec[uid].NNseq >=0)
		NNID=clRec[uid].NNseq;
	else
	  NNID=-1;
	bool taglazy=false;
	bool nrtaglazy=false;
	int id_big=-1;
	REAL errbig=0.0;
	clRec[uid].tagscan=1;
	
	if (id1>=0 && id2>=0)
	{
		if (sq1 >sq2)
		{
			if (clRec[id1].distErr+ clRec[id1].NNdist*sq2/(sq1+sq2) < 0.1*Global::level_step)
			{
				id_big=id1;
				errbig=clRec[id1].distErr+ clRec[id1].NNdist*sq2/(sq1+sq2);
			}
		}
		else
		{
			if (clRec[id2].distErr+ clRec[id2].NNdist*sq1/(sq1+sq2) < 0.1*Global::level_step)
			{
				id_big=id2;
				errbig=clRec[id2].distErr+ clRec[id2].NNdist*sq1/(sq1+sq2);
			}
		}
	}

	//get the nearest leaf of the same branch
	if (thisnode->GetParent()->NumChildren() >1)
	{
		thisnode->GetParent()->ListChildren(tvec);
		int ptr=0;
		
		for (int i=0; i<tvec.size();i++)
		{
			if (tvec[i] !=clRec[uid].Node && tvec[i]->UID !=clRec[uid].NNseq)
			{
				if (QuickDist(id1,id2,tvec[i]->UID,uid,sq1,sq2,dist))
				{
					int NNseq=tvec[i]->UID;
					if (clRec[NNseq].NNdist > dist)
					{
						clRec[NNseq].NNdist=dist;
						if (clRec[NNseq].NNseq >=0)
							clRec[clRec[NNseq].NNseq].NNlist.erase(NNseq);
						clRec[uid].NNlist.insert(NNseq);
						clRec[NNseq].NNseq=uid;
					}
					if (dist < mindist)
					{
						NNID=tvec[i]->UID;
						mindist=dist;
					}
				}
				else
				{
					Global::kdistlist[ptr].idx=i;
					Global::kdistlist[ptr++].dist=clRec[uid].Node->KmerDist(tvec[i]);
				}
			}
		}

		qsort(Global::kdistlist,ptr,sizeof(KDist),CompDist);
	
		for (int i=0;i<ptr;i++)
		{
			if (Global::kdistlist[i].dist < KdistBound(mindist))
			{
				//kdistl.push_back(Global::kdistlist[i].dist);
				if (!QuickDist(id1,id2,tvec[Global::kdistlist[i].idx]->UID,uid,sq1,sq2,dist))
				{
					if (id_big>=0 && Global::cache->Get(id_big,tvec[Global::kdistlist[i].idx]->UID,dist))
					{
						taglazy=true;
						Global::cache->Store(uid,tvec[Global::kdistlist[i].idx]->UID,dist);
					}
					else
					{
						dist=thisnode->NeedleDist(tvec[Global::kdistlist[i].idx]);
						taglazy=false;
					}
				}		
				else
				{
					taglazy=false;
				}
			//	ndistl.push_back(dist);
				
				int NNseq=tvec[Global::kdistlist[i].idx]->UID;

				if (clRec[NNseq].NNdist > dist)
				{
					clRec[NNseq].NNdist=dist;
					if (clRec[NNseq].NNseq >=0)
						clRec[clRec[NNseq].NNseq].NNlist.erase(NNseq);
					clRec[uid].NNlist.insert(NNseq);
					clRec[NNseq].NNseq=uid;
				}
				if (dist < mindist)
				{
					NNID=tvec[Global::kdistlist[i].idx]->UID;
					mindist=dist;
					nrtaglazy=taglazy;
				}
				//mdistl.push_back(mindist);
			}
		}
	}
	
	//enroll all potential NN branches	
	top=thisnode->GetParent();
	
	while (top !=root)
	{
		Tree *tpar=top->GetParent();

		if (top->NumSeqs()> thisnode->NumSeqs() && NNID>=0)
		{
			if (QuickDist(id1,id2,uid,top->UID,sq1,sq2,dist) && dist + mindist < top->GetThres() && !exactsearch) 
			{
				top=tpar;	
				if (!Global::UseApprox)
					continue;
				else
					break;
			}
			if (thisnode->KmerDist(top) < KdistBound(top->GetThres()-mindist)&& !exactsearch)
			{
				top=tpar;	
				if (!Global::UseApprox)
					continue;
				else
					break;
			}
		}	

		if (tpar->NumChildren() >1)
		{
			tvec.clear();	
			tpar->ListChildren(tvec);
			for (int i=0; i<tvec.size();i++)
			{
				if (tvec[i] !=top)
				{
					while (!tvec[i]->IsLeaf() && tvec[i]->NumChildren()==1)
					{
						tvec[i]=tvec[i]->FirstChild();
					}
					dist=thisnode->KmerDist(tvec[i]); 
					if (dist < KdistBound(mindist+tvec[i]->GetThres()))
					{
						heap.Add((void *)tvec[i],dist);
					}
				}
			}
		}
		top=tpar;	
	}
	//perform A* search

	while (!heap.Empty())
	{
		dist=heap.Pop((void *&)top);
 
		if (dist < KdistBound(mindist+top->GetThres()))
		{
			if (top->IsLeaf()) //reach terminal
			{
//				kdistl.push_back(dist);
				if (top->UID==clRec[uid].NNseq) continue;
				if (!QuickDist(id1,id2,top->UID,uid,sq1,sq2,dist))
				{
					if (id_big>=0 && Global::cache->Get(id_big,top->UID,dist))
					{
						taglazy=true;
						Global::cache->Store(top->UID,uid,dist);
					}
					else
					{
						dist=clRec[uid].Node->NeedleDist(top);
						taglazy=false;
					}
				}
				else
				{
					taglazy=false;
				}
				
				int NNseq=top->UID;
				if (clRec[NNseq].NNdist > dist)
				{
					clRec[NNseq].NNdist=dist;
					clRec[uid].NNlist.insert(NNseq);
					if (clRec[NNseq].NNseq >=0)
						clRec[clRec[NNseq].NNseq].NNlist.erase(NNseq);
					clRec[NNseq].NNseq=uid;
				}

				if (dist < mindist)
				{
					NNID=top->UID;
					mindist=dist;
					nrtaglazy=taglazy;
				}
//				mdistl.push_back(mindist);
			}
			else //expand
			{
				tvec.clear();	

				top->ListChildren(tvec);

				for (int i=0; i<tvec.size();i++)
				{
					dist=thisnode->KmerDist(tvec[i]); 
					if (dist < KdistBound(mindist+tvec[i]->GetThres()))
					{
						heap.Add((void *)tvec[i],dist);
					}
				}
			}
		}	
	}
	
/*	if (logon)
	{
		int Numdf=NumAlign1+NumAlign2+NumAlign3-oldal;
		flog << " NumAladded " <<Numdf <<" Seqs " << thisnode->NumSeqs() << endl;
		if (Numdf > 20)
		{
			int sz=mdistl.size();
			for (int i=0;i<sz;i++)
			{
				flog << "(" <<kdistl[i] <<"," <<ndistl[i] <<"," << mdistl[i] <<") ";
			}
			flog <<endl;
		}
	}*/
	clRec[uid].NNdist=mindist;
	if (NNID>=0)
	{
		if (clRec[uid].NNseq >=0 && clRec[uid].NNseq !=NNID)
		{
			clRec[clRec[uid].NNseq].NNlist.erase(uid);
		}
		clRec[uid].NNseq=NNID;
		clRec[NNID].NNlist.insert(uid);

		if (clRec[NNID].NNdist > mindist)
			{
				clRec[NNID].NNdist=mindist;
				clRec[uid].NNlist.insert(NNID);
				if (clRec[NNID].NNseq >=0)
					clRec[clRec[NNID].NNseq].NNlist.erase(NNID);
				clRec[NNID].NNseq=uid;
			}
		if (nrtaglazy)
			clRec[uid].distErr=errbig;
		else
			clRec[uid].distErr=0.0;
	}
}

void ProbCluster::BuildNNlist()
{
	for (int i=0;i<clistTop;i++)
	{
		if (clRec[i].clsid<0)
			NNSearch(i);
	}
}

int MaxHeaps=0;
void ProbCluster::HCluster(fstream &flist,fstream & fstat, fstream &fgroup)
{
	MinHeap clheap;
	MinHeap heapWait;

	REAL dist;
	set<unsigned int>::iterator it;
	set<unsigned int>::iterator it2;	
	
	lastt=clock();
	//logon=true;
	
	REAL DistIsolate=Global::level_max+Global::level_step;
	
/*	logon=true;
	flog.open("debug.log",ios_base::out | ios_base::trunc);	
	Brief(flog);
*/
	for (int i=0;i<clistTop;i++)
		if (clRec[i].clsid < 0)
		{
			if (clRec[i].NNdist <  DistIsolate || Global::ProbSeqOut)
			{	
				clheap.Add((void *)(clRec+i),clRec[i].NNdist);
			}	
			else
			{
				if (clRec[i].NNseq >=0)
					clRec[clRec[i].NNseq].NNlist.erase(i);
				DeleteNode(clRec[i].Node,root);
			}
		}
	
	ClusterRecord *curRec;
	levelcurrent=Global::level_min;
	
	while (!clheap.Empty() || !heapWait.Empty())
	{
		int hsize=clheap.Size();
		if (hsize > MaxHeaps) MaxHeaps=hsize;
		
		
		while (!heapWait.Empty() && (clheap.Empty() || heapWait.Top() < clheap.Top()+Global::level_step/5))
		{
			dist=heapWait.Pop((void *&)curRec);
			//cerr <<" Dequeue " <<curRec->ID << "  dist " << dist <<endl; 
			if (curRec->clsid <0 && curRec->tagscan==0)
			{
				//cerr << "Search NN ... ";
				//if (!Global::UseApprox || curRec->NNseq <0 || clRec[curRec->NNseq].clsid >=0 || (dist > Global::level_min && dist >Global::currentlevel+Global::level_step*0.1))
					NNSearch(curRec->ID);	
			
				//cerr <<" NN " << curRec->NNseq << " dist " << curRec->NNdist <<endl;
			
				if (curRec->NNdist < DistIsolate || Global::ProbSeqOut)
				{
					clheap.Add((void *)curRec,curRec->NNdist);
				}
				else
				{
					if (curRec->NNseq >=0)
					clRec[curRec->NNseq].NNlist.erase(curRec->ID);
			
					DeleteNode(curRec->Node,root);
				}
			}
		}
		
		if (heapWait.Empty() && clheap.Empty())
		{	
			break;
		}
		dist=clheap.Pop((void *&)curRec);
	
		if (curRec->clsid>= 0 || curRec->NNseq<0 || clRec[curRec->NNseq].clsid >= 0) continue;
		if (curRec->NNdist > dist) continue; //outdated
			


		while (dist >levelcurrent) 
		{
			Report(flist,fstat,fgroup);
			levelcurrent+=Global::level_step;
			if (levelcurrent > Global::level_max+Global::level_step*0.1)
			{
			//	if (logon) flog.close();
				return;	
			}
		}
		
		if (NumCls<=1) 
		{
			Report(flist,fstat,fgroup);
		//	if (logon) flog.close();
				return;	
		}
		
		Global::currentlevel=dist;
	
		int rNN=clRec[curRec->NNseq].NNseq;
		if (rNN>=0 && rNN !=curRec->ID)
			clRec[rNN].NNlist.erase(curRec->NNseq);
		//cerr << curRec->ID << " + " << curRec->NNseq << "-->" << clistTop <<endl;

		Tree *newnode=MergeNode(curRec->Node,clRec[curRec->NNseq].Node,clistTop);
		Tree *top=FindCommonParent(curRec->Node,clRec[curRec->NNseq].Node);

		DeleteNode(curRec->Node,top);
		DeleteNode(clRec[curRec->NNseq].Node,top);
		
		AddRecord(clistTop,newnode);

		curRec->NNlist.erase(curRec->NNseq);
		clRec[curRec->NNseq].NNlist.erase(curRec->ID);
		
		for (it=curRec->NNlist.begin() ; it != curRec->NNlist.end(); it++ )
		{
			REAL cdist,rdist;
			REAL mdist=maxdist;
			int minid=-1;
			
			clRec[*it].NNlist.erase(curRec->ID);
			clRec[*it].NNlist.erase(curRec->NNseq);
			
			heapWait.Add((void *)(clRec+(*it)),clRec[*it].NNdist);

			for (it2=clRec[*it].NNlist.begin();it2!=clRec[*it].NNlist.end();it2++)
			{
				if (clRec[*it2].NNdist <mdist)
				{
					mdist=clRec[*it2].NNdist;
					minid=*it2;
				}
			}
			   
			if (Global::cache->Get(clRec[*it].ID,curRec->NNseq,cdist))
			{
				if (Global::Dist_Style == DS_Ave)
				{
					rdist=(curRec->NumSeqs*clRec[*it].NNdist+clRec[curRec->NNseq].NumSeqs*cdist)/(curRec->NumSeqs+clRec[curRec->NNseq].NumSeqs);
				}
				else
				{
					rdist=max(clRec[*it].NNdist,cdist);
				}
				
				Global::cache->Store(clistTop,clRec[*it].ID,rdist);
				if (rdist <mdist)
				{
					mdist=rdist;
					minid=clistTop;
				}	
				if (clRec[clistTop].NNdist > rdist)
				{
					clRec[clistTop].NNseq=*it;
					clRec[clistTop].NNdist = rdist;
				}
			}
			
			clRec[*it].NNdist=mdist;
			clRec[*it].NNseq=minid;
			clRec[*it].tagscan=0;
			
			if (minid >=0)
				clRec[minid].NNlist.insert(*it);
			
		}
		
		for (it=clRec[curRec->NNseq].NNlist.begin() ; it != clRec[curRec->NNseq].NNlist.end(); it++ )
		{
			REAL cdist,rdist;
			REAL mdist=maxdist;
			int minid=-1;
			
			clRec[*it].NNlist.erase(curRec->ID);
			clRec[*it].NNlist.erase(curRec->NNseq);

			heapWait.Add((void *)(clRec+(*it)),clRec[*it].NNdist);
		
			for (it2=clRec[*it].NNlist.begin();it2!=clRec[*it].NNlist.end();it2++)
			{
				if (clRec[*it2].NNdist <mdist)
				{
					mdist=clRec[*it2].NNdist;
					minid=*it2;
				}
			}

			if (Global::cache->Get(clRec[*it].ID,curRec->ID,cdist))
			{
				if (Global::Dist_Style == DS_Ave)
				{
					rdist=(curRec->NumSeqs*cdist+clRec[curRec->NNseq].NumSeqs*clRec[*it].NNdist)/(curRec->NumSeqs+clRec[curRec->NNseq].NumSeqs);
				}
				else
				{
					rdist=max(clRec[*it].NNdist,cdist);
				}
				
				Global::cache->Store(clistTop,clRec[*it].ID,rdist);
				
				if (rdist < mdist)
				{
					mdist=rdist;
					minid=clistTop;
				}	
				if (clRec[clistTop].NNdist > rdist)
				{
					clRec[clistTop].NNseq=*it;
					clRec[clistTop].NNdist = rdist;
				}
			}

			clRec[*it].NNdist=mdist;
			clRec[*it].NNseq=minid;
			clRec[*it].tagscan=0;

			if (minid >=0)
				clRec[minid].NNlist.insert(*it);
		}
		
		if (clRec[clistTop].NNseq >=0)
		{	
				clRec[clRec[clistTop].NNseq].NNlist.insert(clistTop);
		}		

		curRec->NNlist.clear();
		clRec[curRec->NNseq].NNlist.clear();
		curRec->clsid=clistTop;
		clRec[curRec->NNseq].clsid=clistTop;


		clRec[clistTop].Seqlist=curRec->Seqlist;
		clRec[clistTop].Seqlist.insert(clRec[curRec->NNseq].Seqlist.begin(),clRec[curRec->NNseq].Seqlist.end());
		curRec->Seqlist.clear();
		clRec[curRec->NNseq].Seqlist.clear();
		
		AddNode(newnode,top);
		

		NNSearch(clistTop,true,curRec->ID,curRec->NNseq,curRec->NumSeqs,clRec[curRec->NNseq].NumSeqs);
		
		if (clRec[clistTop].NNdist < DistIsolate || Global::ProbSeqOut)
		{
			clheap.Add((void *)(clRec+clistTop),clRec[clistTop].NNdist);
		}
		else
		{
			if (clRec[clistTop].NNseq >=0)
				clRec[clRec[clistTop].NNseq].NNlist.erase(clistTop);
			DeleteNode(clRec[clistTop].Node,root);
		}
		clistTop++;
		NumCls--;
	}
	
	while (levelcurrent < Global::level_max) 
		{
			Report(flist,fstat,fgroup);
			levelcurrent+=Global::level_step;
		}
		
}

void ProbCluster::Report(fstream &flist,fstream &fstat,fstream &fgroup)
{
	int i,top;
	set <unsigned int>::iterator it;
	fstream fclust;

	thist=clock();	
	printf("%6.3f %d    Als %d Kmer %ld Elapse %10.3g s\n",levelcurrent,NumCls,NumAlign1+NumAlign2+NumAlign3,
	Kcomp+Kcomp2,(REAL)(thist-lastt)/CLOCKS_PER_SEC);
	
	fgroup.precision(3);
	flist.precision(3);
	fstat.precision(3);
	lastt=thist;
	
	if(Global::ProbSeqOut)
	{
		char buf[2000];
		sprintf(buf,"%s.%3g.psq",Global::ProbSeqOut,levelcurrent);
		fclust.open(buf,ios_base::out | ios_base::trunc);	
	}
	
	flist << fixed << levelcurrent << " " << NumCls << endl;
	fgroup << fixed << levelcurrent << " |";
	fstat << fixed << levelcurrent << " ";
	
	for (i=0;i<clistTop;i++)
		{
			if (clRec[i].clsid<0)
			{
				fstat << fixed << clRec[i].NumSeqs<<" ";
				for (it=clRec[i].Seqlist.begin();it!=clRec[i].Seqlist.end();it++)
				{
					if (it !=clRec[i].Seqlist.begin())
						fgroup <<" ";
					fgroup << SeqID[*it];
				}
				fgroup << "|";
				if (Global::ProbSeqOut)
				{
					if (clRec[i].Node->SingleSeq)
					{
						ProbString pr(SeqStrs[clRec[i].ID]);
						//fclust << clRec[i].ID << endl;
						fclust << pr;
					}
					else
					{
						//fclust << clRec[i].ID << endl;
						fclust << clRec[i].Node->Center();
					}
				}

			}
		}

	fstat << endl;
	fgroup <<endl;	
	if (Global::ProbSeqOut)
		fclust.close();
}

void ProbCluster::PrintNN(fstream & out, int SeqNum)
{ 
	set<unsigned int>::iterator it;
	for (int i=0;i<SeqNum;i++)
	{ 	
		out << " ID " << i << " NN " << clRec[i].NNseq << " dist " <<clRec[i].NNdist << " clsid "<< clRec[i].clsid;
		out << " NNlist ";
		for (it=clRec[i].NNlist.begin();it!=clRec[i].NNlist.end();it++)
			out << *it <<" ";
		out << endl;
	}
}


/*
void ProbCluster::ReportList(fstream &flist,fstream & fstat, fstream &fgroup, REAL leveldown, REAL levelup)
{
	vector<Tree *> repvec;
	vector<Tree *> leafvec;
	REAL level=leveldown;
	int gsize,lsize,i,j,numseq;
	Tree *current,*leaf;
	
	while (level <=levelup*1.001)
	{
		root->ListChildrenAt(level, repvec);
		gsize=repvec.size();
		fgroup.precision(2);
		flist.precision(2);
		fstat.precision(2);
		flist << fixed << level << " " << gsize << endl;	
		fgroup << fixed << level << " ";
		fstat << fixed << level << " ";
		
		for (i=0;i<gsize;i++)
		{
			numseq=0;
			current=repvec[i];
			current->ListLeaf(leafvec);
			lsize=leafvec.size();
			for (j=0;j<lsize;j++)
			{
				leaf=leafvec[j];
				fgroup << leaf->Center().UId();
				numseq+=leaf->NumSeqs();
				if (j < lsize-1)
					fgroup << ' ';
				else
					fgroup <<'|';
			}
			leafvec.clear();
			fstat << numseq << " ";
		}
		fstat << endl;
		fgroup << endl;
		repvec.clear();
		level*=levelinc;
	}
}

*/
