#include "stdinc.h"
#include "newhash.h"
#include "extfunc.h"
#include "extvab.h"
#include "dfibHeap.h"
#include "fibHeap.h"
#include "darray.h"


#define CNBLOCKSIZE 10000
#define GAPARRSIZE 256
#define BIG_NEG -10000000
#define BIG_POS 10000000
static STACK *isStack;
static int onsameCtgPE;
extern int calcuIS(STACK *intStack, int *SD);
void outputBundle(FILE *fp, int insertS);

static CONNECT *bun1AccuConnect(unsigned int e1, unsigned int e2, int gap, int weight)
{
  if(e1 == e2 || e1 == getTwinCtg(e2))
    return NULL;

  CONNECT *connect = NULL;

  //long long sum;
  if(weight > 255)
    weight = 255;

  connect = getCntBetween(e1, e2);

  if(connect)
    {
      if(!weight)
        return connect;

      //existCounter++;
      //if(!inherit){
      //sum = connect->weightNotInherit*connect->gapLen + gap*weight;
      //connect->gapLen = sum/(connect->weightNotInherit+weight);
      int i = connect->weightNotInherit;

      if(connect->weightNotInherit + weight <= 255)
        connect->weightNotInherit += weight;
      else if(connect->weightNotInherit < 255)
        connect->weightNotInherit = 255;

      for(; i < connect->weightNotInherit; i++)
        {
          //connect->PE[i]=gap;
          //fprintf(stderr,"inputting a PE with estimated gap size %d\n",gap);
        }

      /*}else{
      	//sum = connect->weight*connect->gapLen + gap*weight;
      	//connect->gapLen = sum/(connect->weight+weight);
      	if(!connect->inherit){
      		connect->maxSingleWeight = connect->weightNotInherit;
      	}
      	connect->inherit = 1;
      	connect->maxSingleWeight = connect->maxSingleWeight>weight ?
      			connect->maxSingleWeight:weight;
      }*/
      if(connect->weight + weight <= 255)
        {
          connect->weight += weight;
        }
      else if(connect->weight < 255)
        {
          connect->weight = 255;
        }

    }
  else
    {
      newCntCounter++;
      connect = allocateCN(e2, gap);

      if(cntLookupTable)
        putCnt2LookupTable(e1, connect);

      connect->weight = weight;
      //connect->PE=(int *)ckalloc(GAPARRSIZE*sizeof(int));//newly added
      //fprintf(stderr,"creating array for PEs in a connection.\n");
      int i;

      for(i = 0; i < weight; i++)
        {
          //connect->PE[i]=gap;
          //fprintf(stderr,"inputting a PE with estimated gap size %d\n",gap);
        }

      if(contig_array[e1].mask || contig_array[e2].mask)
        {
          connect->mask = 1;
        }

      connect->next = contig_array[e1].downwardConnect;
      contig_array[e1].downwardConnect = connect;
      //if(!inherit){
      connect->weightNotInherit = weight;
      /*}else{
      	connect->weightNotInherit = 0;
      	connect->inherit = 1;
      	connect->maxSingleWeight = weight;
      }*/
    }

  return connect;
}

static int in1PE(unsigned int e1, int pre_pos, unsigned int bal_e2, int pos, int insert_size)
{
  int gap, realpeSize;
  unsigned int bal_e1, e2;

  if(e1 == bal_e2)
    {
      //ignorePE1++;
      return -1;  //orientation wrong
    }

  bal_e1 = getTwinCtg(e1);
  e2 = getTwinCtg(bal_e2);

  if(e1 == e2)
    {
      realpeSize = contig_array[e1].length + overlaplen - pre_pos - pos;

      if(realpeSize > 0)
        {
          //peSUM += realpeSize;
          onsameCtgPE++;

          if((int)contig_array[e1].length > insert_size)
            {
              int *item = (int *)stackPush(isStack);
              (*item) = realpeSize;
            }
        }

      return 2;
    }

  gap = insert_size - overlaplen + pre_pos + pos - contig_array[e1].length - contig_array[e2].length;
  //fprintf(stderr,"[%s]\t%d\t%d\tgap\t%d\t%d\t%d\t%d\n",__FUNCTION__,e1,e2,gap,contig_array[e1].bal_edge,contig_array[e2].bal_edge,insert_size);
  //if(gap<-(insert_size/10)){
  //	//ignorePE2++;
  //	return 0;
  //}
  bun1AccuConnect(e1, e2, gap, 1);
  bun1AccuConnect(bal_e2, bal_e1, gap, 1);

  return 1;
}

static int inputPE(FILE *fp, int peGrad, char *line)
{
  long long pre_readno, readno, minno, maxno;
  int pre_pos, pos, flag, PE, count = 0;
  unsigned int pre_contigno, contigno, newIndex;

  if(peGrad < 0 || peGrad > gradsCounter)
    {
      printf("[%s]specified pe grad is out of bound .\n", __FUNCTION__);
      return 0;
    }

  maxno = pes[peGrad].PE_bound;

  if(peGrad == 0)
    minno = 0;
  else
    minno = pes[peGrad - 1].PE_bound;

  //onsameCtgPE = peSUM = 0;
  PE = pes[peGrad].insertS;

  if(strlen(line))
    {
      sscanf(line, "%lld %d %d", &pre_readno, &pre_contigno, &pre_pos);

      //printf("first record %d %d %d\n",pre_readno,pre_contigno,pre_pos);
      if(pre_readno <= minno)
        pre_readno = -1;
    }
  else
    pre_readno = -1;

  //ignorePE1 = ignorePE2 = ignorePE3 = ignorePE4 = ignorePE5 = 0;
  //static_flag = 1;
  isStack = (STACK *)createStack(CNBLOCKSIZE, sizeof(int));

  while(fgets(line, lineLen, fp) != NULL)
    {
      sscanf(line, "%lld %d %d", &readno, &contigno, &pos);

      if(readno > maxno)
        break;

      if(readno <= minno)
        continue;

      newIndex = index_array[contigno];

      //if(contig_array[newIndex].bal_edge==0)
      if(isSameAsTwin(newIndex))
        continue;

      if(PE && (readno % 2 == 0) && (pre_readno == readno - 1)) // they are a pair of reads
        {
          flag = in1PE(pre_contigno, pre_pos, newIndex, pos, PE);

          if(flag == 1)
            count++;
        }

      pre_readno = readno;
      pre_contigno = newIndex;
      pre_pos = pos;
    }

  printf("[%s]Finish loading all PEs in grad %d .\n", __FUNCTION__, peGrad);
  printf("[%s]Calculating estimated gap size for all connections .\n", __FUNCTION__);
  /*unsigned int i;
  for(i=1;i<=num_ctg;i++){
  		CONNECT *tmp=contig_array[i].downwardConnect;
  		while(tmp){
  			if(tmp->weightNotInherit<=8&&tmp->weightNotInherit>2){//delete max and min value
  				int max=BIG_NEG,maxid=-1,min=BIG_POS,minid=-1;
  				int weight=tmp->weightNotInherit;
  				int ii;
  				for(ii=0;ii<weight;ii++){
  					if(tmp->PE[ii]>max){
  						max=tmp->PE[ii];
  						maxid=ii;
  					}
  					if(tmp->PE[ii]<=min){
  						min=tmp->PE[ii];
  						minid=ii;
  					}
  				}
  				int sum=0;
  				for(ii=0;ii<weight;ii++){
  					if(ii!=maxid&&ii!=minid){
  						sum+=tmp->PE[ii];
  					}
  				}
  				//ignorePE4+=2;
  				tmp->gapLen=sum/(weight-2);
  				//fprintf(stderr,"estimating contigs' gap by removing max&min PE ,with max&min %d %d\n",
  					//tmp->PE[maxid],tmp->PE[minid]);
  			}else if(tmp->weightNotInherit>8){//delete values exceed 3*SD
  				long long int sum=0;
  				int weight=tmp->weightNotInherit;
  				int ii;
  				int counter=0;
  				for(ii=0;ii<weight;ii++){
  					sum+=tmp->PE[ii];
  				}

  				long long int avg=sum/weight;
  				sum = 0;
  				for(ii=0;ii<weight;ii++){
  					sum+=((avg-(long long int)tmp->PE[ii])*(avg-(long long int)tmp->PE[ii]));
  				}

  				double SD=(sqrt((double)sum/(weight-1)))*3;//just for fast
  				sum=0;
  				int num=0;
  				for(ii=0;ii<weight;ii++){
  					if(abs(tmp->PE[ii]-avg)<=SD){
  						sum+=tmp->PE[ii];
  						num++;
  					}else{
  						//ignorePE5++;
  						counter++;
  					}
  				}
  				if(num==0){
  					//fprintf(stderr,"[%s]num=0 in removing exceed 3*SD(%.1f) avg(%lld)step",__FUNCTION__,SD,avg);
  					for(ii=0;ii<weight;ii++){
  						fprintf(stderr,"%d\t",tmp->PE[ii]);
  					}
  				}
  				tmp->gapLen=sum/num;
  				//fprintf(stderr,"estimating contigs' gap by removing PE exceeding 3*SD ,removing %d PEs\n",counter);
  			}else if(tmp->weightNotInherit<=2){
  				int weight=tmp->weightNotInherit;
  				int sum=0;
  				int ii;
  				for(ii=0;ii<weight;ii++){
  					sum+=tmp->PE[ii];
  				}
  				tmp->gapLen=sum/weight;
  				//fprintf(stderr,"weight too small , directly estimate gap size.\n");
  			//}
  			//fprintf(stderr,"finish %d connection.\n",i);
  			//free((void *)tmp->PE);
  			tmp=tmp->next;
  		}
  }*/
  //printf("%d PEs with insert size %d attached, %d + %d + %d ignored\n",count,PE,ignorePE1,ignorePE2,ignorePE3);
  fprintf(stderr, "[%s]%d PEs of insert size %d loaded .\n", __FUNCTION__, count, PE);
  //fprintf(stderr,"[%s]PEs discarded:%d because of wrong orientation,%d too close,%d too far,\n",__FUNCTION__,ignorePE1,ignorePE2,ignorePE3);
  //fprintf(stderr,"[%s]%d deleted by removing max&min , %d not fall in 3*SD.\n",__FUNCTION__,ignorePE4,ignorePE5);
  //printf("[%s]%d PEs of insert size %d loaded .\n",__FUNCTION__,count,PE);
  //printf("[%s]PEs discarded :%d because of wrong orientation,%d too close,%d too far ,\n",__FUNCTION__,ignorePE1,ignorePE2,ignorePE3);
  //printf("[%s]%d deleted by removing max&min , %d not fall in 3*SD .\n",__FUNCTION__,ignorePE4,ignorePE5);

  /*if(onsameCtgPE>0){
  	//printf("estimated PE size %lli, by %d pairs\n",peSUM/onsameCtgPE,onsameCtgPE);
  	int SD=0;
  	int avg=calcuIS(isStack,&SD);
  	printf("[%s]%d PE attached on same contig with estimated insert size %d SD %d .\n",__FUNCTION__,onsameCtgPE,avg,SD);
  }*/
  //printf("on contigs longer than %d, %d pairs found,",PE,isStack->item_c);
  //printf("insert_size estimated: %d\n",calcuIS(isStack));
  //freeStack(isStack);
  return count;
}

int call_bundle()
{
  char name[256], *line;
  FILE *fp, *linkF;
  int i;
  int flag = 0;
  unsigned int j;

  loadUpdatedEdges(graphfile);

  //sprintf(name,"%s.bundle",graphfile);

  linkF = ckopen(name, "w");

  if(!pes)
    loadPEgrads(graphfile);

  sprintf(name, "%s.readOnContig", graphfile);
  fp = ckopen(name, "r");

  lineLen = 1024;
  line = (char *)ckalloc(lineLen * sizeof(char));

  fgets(line, lineLen, fp);
  line[0] = '\0';

  //printf("\n");
  newCntCounter = 0;

  //createCntMemManager();
  //createCntLookupTable();
  /*int *length_array = (unsigned int *)ckalloc((num_ctg+1)*sizeof(unsigned int));
  //use length_array to change info in index_array
  for(i=1;i<=num_ctg;i++)
  	length_array[i] = 0;

  for(i=1;i<=num_ctg;i++){
  	if(index_array[i]>0)
  		length_array[index_array[i]] = i;
  }
  for(i=1;i<=num_ctg;i++)
  	index_array[i] = length_array[i];
  */
  for(i = 0; i < gradsCounter; i++)
    {

      createCntMemManager();
      createCntLookupTable();
      //
      flag += inputPE(fp, i, line);
      //sprintf(name,"%d.bundle",i);

      //printf("%lld new connections\n",newCntCounter/2);
      /*if(!flag){
      	destroyConnectMem();
      	deleteCntLookupTable();
      	for(j=1;j<=num_ctg;j++)
      		contig_array[j].downwardConnect = NULL;
      	//printf("\n");
      	continue;
      }*/
      flag = 0;
      //linkF= ckopen(name,"w");
      //outputBundle(linkF, pes[i].insertS);

      for(j = 1; j <= num_ctg; j++)
        {
          CONNECT *tmp = contig_array[j].downwardConnect;

          while(tmp)
            {
              free((void *)tmp->PE);
              tmp = tmp->next;
            }

          contig_array[j].downwardConnect = NULL;
        }

      //destroyConnectMem();
      //deleteCntLookupTable();

      fclose(linkF);
    }

  outputBundle(linkF, 1);
  destroyConnectMem();
  deleteCntLookupTable();

  free((void *)line);
  fclose(fp);
  //fclose(linkF);
  printf("[%s]all PEs attached\n", __FUNCTION__);

  return 0;
}

void outputBundle(FILE *fp, int insertS)
{
  unsigned int i, bal_ctg, bal_toCtg;
  CONNECT *cnts, *temp_cnt;

  //printf("outputLinks, %d contigs\n",num_ctg);
  for(i = 1; i <= num_ctg; i++)
    {
      cnts = contig_array[i].downwardConnect;
      bal_ctg = getTwinCtg(i);

      //fprintf(stderr,"contig %d.\n",i);
      while(cnts)
        {
          if(cnts->weightNotInherit <= bund_threshold)
            {
              cnts = cnts->next;
              continue;
            }

          //fprintf(stderr,"with contig %d.\n",cnts->contigID);
          //fprintf(fp,"%-10d %-10d\t%d\t%d\t%d\n"
          //,i,cnts->contigID,cnts->gapLen,cnts->weight,insertS);
          /*int st1,st2,ed1,ed2,len1,len2,gap;
          len1=contig_array[i].length+overlaplen;
          len2=contig_array[cnts->contigID].length+overlaplen;
          gap=-cnts->gapLen;
          if(len1<gap){
          	st1=0;
          	if(len2<gap){
          		ed1=len1-gap+len2;
          		st2=gap-len1;
          		ed2=len2;
          	}else{
          		ed1=len1;
          		st2=gap-len1;
          		ed2=gap;
          }
          }else{
          	st1=len1+overlaplen-gap;
          	st2=0;
          	if(len2<gap){
          		ed1=st1+len2;
          		ed2=len2;
          	}else{
          		ed2=gap;
          		ed1=len1;
          	}
          }

          unsigned int id1,id2;
          id1=index_array[i];
          id2=index_array[cnts->contigID];*/
          /*if((id1/2+1)==1194){
          	int ii;
          	fprintf(stdout,"\n");
          	for(ii=0;ii<cnts->weightNotInherit;++ii){
          		fprintf(stdout,"%d ",cnts->PE[ii]);
          	}
          	fprintf(stdout,"\n");
          }*/
          /*if(isSmallerThanTwin(id1)){
          	if(isSmallerThanTwin(id2)){
          		fprintf(fp,"%u\t%d\t%u\t%d\t%d\n",id1/2+1,len1,id2/2+1,len2,cnts->gapLen,cnts->weightNotInherit);

          	}else{
          		fprintf(fp,"%u\t%d\t%u\t%d\t%d\n",id1/2+1,len1,-id2/2,len2,cnts->gapLen,cnts->weightNotInherit);
          	}
          }else{
          	if(isSmallerThanTwin(id2)){
          		fprintf(fp,"%u\t%d\t%u\t%d\t%d\n",-id1/2,len1,id2/2+1,len2,cnts->gapLen,cnts->weightNotInherit);
          	}else{
          		fprintf(fp,"%u\t%d\t%u\t%d\t%d\n",-id1/2,len1,-id2/2,len2,cnts->gapLen,cnts->weightNotInherit);
          	}
          }*/
          //int ii=0;
          //int weight=cnts->weightNotInherit;
          //for(;ii<weight;++ii){
          //	fprintf(fp,"%d\t%d\t%d\t",icnts->gapLen);
          //}
          if(cnts->gapLen < 0)
            {
              fprintf(fp, "%d\t%d\t%d\n", i, cnts->contigID, cnts->gapLen);
            }

          //fprintf(fp,"\n");
          cnts->weightNotInherit = 0;

          bal_toCtg = getTwinCtg(cnts->contigID);
          temp_cnt = getCntBetween(bal_toCtg, bal_ctg);

          if(temp_cnt)
            temp_cnt->weightNotInherit = 0;

          cnts = cnts->next;
        }
    }
}