#include "stdinc.h"
#include "newhash.h"
#include "extfunc.h"
#include "extvab.h"

#define RDBLOCKSIZE 50
#define CTGappend 50

static Kmer MAXKMER;
static int Ncounter;
static int allGaps;

// for multi threads
static int *counters;
static pthread_mutex_t mutex;
static int scafBufSize = 100;
static boolean *flagBuf;
static unsigned char *thrdNoBuf;
static STACK **ctgStackBuffer;
static int scafCounter;
static int scafInBuf;

static void MarkCtgOccu(unsigned int ctg);
/*
static void printRead(int len,char *seq)
{
	int j;
	fprintf(stderr,">read\n");
	for(j=0;j<len;j++)
		fprintf(stderr,"%c",int2base((int)getCharInTightString(seq,j)));
	fprintf(stderr,"\n");
}
*/
static void attach1read2contig(unsigned int ctgID, int len, int pos, long long starter)
{
  unsigned int ctg = index_array[ctgID]; //new index in contig array

  if(isLargerThanTwin(ctg))
    {
      ctg = getTwinCtg(ctg);  // put all reads in one contig of a twin
      pos = contig_array[ctg].length + overlaplen - pos - len;
    }

  if(!contig_array[ctg].closeReads)
    contig_array[ctg].closeReads = (STACK *)createStack(RDBLOCKSIZE, sizeof(READNEARBY));

  READNEARBY *rd = (READNEARBY *)stackPush(contig_array[ctg].closeReads);
  rd->len = len;
  rd->dis = pos;
  rd->seqStarter = starter;
}

static void convertIndex()
{
  int *length_array = (int *)ckalloc((num_ctg + 1) * sizeof(int));
  unsigned int i;

  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];  //contig i with new index: index_array[i]

  free((void *)length_array);

}

static long long getRead1by1(FILE *fp, DARRAY *readSeqInGap)
{
  long long readCounter = 0;

  if(!fp)
    return readCounter;

  int len, ctgID, pos;
  long long starter;
  char *pt;
  char *freadBuf = (char *)ckalloc((maxReadLen / 4 + 1) * sizeof(char));

  while(fread(&len, sizeof(int), 1, fp) == 1)
    {
      if(fread(&ctgID, sizeof(int), 1, fp) != 1)
        break;

      if(fread(&pos, sizeof(int), 1, fp) != 1)
        break;

      if(fread(freadBuf, sizeof(char), len / 4 + 1, fp) != (unsigned)(len / 4 + 1))
        break;

      //put seq to dynamic array
      starter = readSeqInGap->item_c;

      if(!darrayPut(readSeqInGap, starter + len / 4)) // make sure there's room for this seq
        break;

      pt = (char *)darrayPut(readSeqInGap, starter);
      bcopy(freadBuf, pt, len / 4 + 1);
      attach1read2contig(ctgID, len, pos, starter);
      readCounter++;
    }

  free((void *)freadBuf);
  return readCounter;
}
// Darray *readSeqInGap
static boolean loadReads4gap(char *graphfile)
{
  FILE *fp, *fp2;
  char name[1024];
  long long readCounter;

  sprintf(name, "%s.readInGap", graphfile);
  fp = fopen(name, "rb");
  sprintf(name, "%s.longReadInGap", graphfile);
  fp2 = fopen(name, "rb");

  if(!fp && !fp2)
    return 0;

  if(!orig2new)
    {
      convertIndex();
      orig2new = 1;
    }

  readSeqInGap = (DARRAY *)createDarray(1000000, sizeof(char));

  if(fp)
    {
      readCounter = getRead1by1(fp, readSeqInGap);
      //printf("Loaded %lld reads from %s.readInGap\n",readCounter,graphfile);
      fclose(fp);
    }

  if(fp2)
    {
      readCounter = getRead1by1(fp2, readSeqInGap);
      //printf("Loaded %lld reads from %s.LongReadInGap\n",readCounter,graphfile);
      fclose(fp2);
    }

  return 1;
}

static void debugging1()
{
  unsigned int i;

  if(orig2new)
    {
      unsigned 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];  //contig i with original index: index_array[i]

      orig2new = 0;
    }

  READNEARBY *rd;
  int j;
  char *pt;

  for(i = 1; i <= num_ctg; i++)
    {
      if(!contig_array[i].closeReads)
        continue;

      if(index_array[i] != 735)
        continue;

      //printf("contig %d, len %d: \n",index_array[i],contig_array[i].length);
      stackBackup(contig_array[i].closeReads);

      while((rd = (READNEARBY *)stackPop(contig_array[i].closeReads)) != NULL)
        {
          printf("%d\t%d\t%lld\t", rd->dis, rd->len, rd->seqStarter);
          pt = (char *)darrayGet(readSeqInGap, rd->seqStarter);

          for(j = 0; j < rd->len; j++)
            printf("%c", int2base((int)getCharInTightString(pt, j)));

          printf("\n");
        }

      stackRecover(contig_array[i].closeReads);
    }

}

static void initiateCtgInScaf(CTGinSCAF *actg)
{
  actg->cutTail = 0;
  actg->cutHead = overlaplen;
  actg->gapSeqLen = 0;
}

static int procGap(char *line, STACK *ctgsStack)
{
  char *tp;
  int length, i, seg;
  unsigned int ctg;
  CTGinSCAF *ctgPt;

  tp = strtok(line, " ");
  tp = strtok(NULL, " "); //length
  length = atoi(tp);
  tp = strtok(NULL, " "); //seg
  seg = atoi(tp);

  if(!seg)
    return length;

  for(i = 0; i < seg; i++)
    {
      tp = strtok(NULL, " ");
      ctg = atoi(tp);
      MarkCtgOccu(ctg);
      ctgPt = (CTGinSCAF *)stackPush(ctgsStack);
      initiateCtgInScaf(ctgPt);
      ctgPt->ctgID = ctg;
      ctgPt->start = 0;
      ctgPt->end = 0;
      ctgPt->scaftig_start = 0;
      ctgPt->mask = 1;
    }

  return length;
}

static void debugging2(int index, STACK *ctgsStack)
{
  CTGinSCAF *actg;

  stackBackup(ctgsStack);
  printf(">scaffold%d\t%d 0.0\n", index, ctgsStack->item_c);

  while((actg = stackPop(ctgsStack)) != NULL)
    {
      printf("%d\t%d\t%d\t%d\n",
             actg->ctgID, actg->start, actg->end, actg->scaftig_start);
    }

  stackRecover(ctgsStack);
}

static int cmp_reads(const void *a, const void *b)
{
  READNEARBY *A, *B;
  A = (READNEARBY *)a;
  B = (READNEARBY *)b;

  if(A->dis > B->dis)
    return 1;
  else if(A->dis == B->dis)
    return 0;
  else
    return -1;
}

static void cutRdArray(READNEARBY *rdArray, int gapStart, int gapEnd, int *count, int arrayLen, READNEARBY *cutArray)
{
  int i;
  int num = 0;

  for(i = 0; i < arrayLen; i++)
    {
      if(rdArray[i].dis > gapEnd)
        break;

      if((rdArray[i].dis + rdArray[i].len) >= gapStart)
        {
          cutArray[num].dis = rdArray[i].dis;
          cutArray[num].len = rdArray[i].len;
          cutArray[num++].seqStarter = rdArray[i].seqStarter;
        }
    }

  *count = num;
}

static void outputTightStr(FILE *fp, char *tightStr, int start, int length, int outputlen, int revS, int *col)
{
  int i;
  int end;
  int column = *col;

  if(!revS)
    {
      end = start + outputlen <= length ? start + outputlen : length;

      for(i = start; i < end; i++)
        {
          fprintf(fp, "%c", int2base((int)getCharInTightString(tightStr, i)));

          if((++column) % 100 == 0)
            {
              //column = 0;
              fprintf(fp, "\n");
            }
        }
    }
  else
    {
      end = length - start - outputlen - 1 >= 0 ? length - start - outputlen : 0;

      for(i = length - 1 - start; i >= end; i--)
        {
          fprintf(fp, "%c", int2compbase((int)getCharInTightString(tightStr, i)));

          if((++column) % 100 == 0)
            {
              fprintf(fp, "\n");
              //column = 0;
            }
        }
    }

  *col = column;
}

static void outputTightStrLowerCase(FILE *fp, char *tightStr, int start, int length, int outputlen, int revS, int *col)
{
  int i;
  int end;
  int column = *col;

  if(!revS)
    {
      end = start + outputlen <= length ? start + outputlen : length;

      for(i = start; i < end; i++)
        {
          fprintf(fp, "%c", "actg"[(int)getCharInTightString(tightStr, i)]);

          if((++column) % 100 == 0)
            {
              //column = 0;
              fprintf(fp, "\n");
            }
        }
    }
  else
    {
      end = length - start - outputlen - 1 >= 0 ? length - start - outputlen : 0;

      for(i = length - 1 - start; i >= end; i--)
        {
          fprintf(fp, "%c", "tgac"[(int)getCharInTightString(tightStr, i)]);

          if((++column) % 100 == 0)
            {
              fprintf(fp, "\n");
              //column = 0;
            }
        }
    }

  *col = column;
}

static void outputNs(FILE *fp, int gapN, int *col)
{
  int i, column = *col;

  for(i = 0; i < gapN; i++)
    {
      fprintf(fp, "N");

      if((++column) % 100 == 0)
        {
          //column = 0;
          fprintf(fp, "\n");
        }
    }

  *col = column;
}

static void outputGapInfo(unsigned int ctg1, unsigned int ctg2)
{
  unsigned int bal_ctg1 = getTwinCtg(ctg1);
  unsigned int bal_ctg2 = getTwinCtg(ctg2);

  if(isLargerThanTwin(ctg1))
    fprintf(stderr, "%d\t", index_array[bal_ctg1]);
  else
    fprintf(stderr, "%d\t", index_array[ctg1]);

  if(isLargerThanTwin(ctg2))
    fprintf(stderr, "%d\n", index_array[bal_ctg2]);
  else
    fprintf(stderr, "%d\n", index_array[ctg2]);
}

static void output1gap(FILE *fo, int scafIndex, CTGinSCAF *prevCtg, CTGinSCAF *actg, DARRAY *gapSeqArray)
{
  unsigned int ctg1, bal_ctg1, length1;
  int start1, outputlen1;
  unsigned int ctg2, bal_ctg2, length2;
  int start2, outputlen2;
  char *pt;
  int column = 0;

  ctg1 = prevCtg->ctgID;
  bal_ctg1 = getTwinCtg(ctg1);
  start1 = prevCtg->cutHead;
  length1 = contig_array[ctg1].length + overlaplen;

  if(length1 - prevCtg->cutTail - start1 > CTGappend)
    {
      outputlen1 = CTGappend;
      start1 = length1 - prevCtg->cutTail - outputlen1;
    }
  else
    outputlen1 = length1 - prevCtg->cutTail - start1;

  ctg2 = actg->ctgID;
  bal_ctg2 = getTwinCtg(ctg2);
  start2 = actg->cutHead;
  length2 = contig_array[ctg2].length + overlaplen;

  if(length2 - actg->cutTail - start2 > CTGappend)
    {
      outputlen2 = CTGappend;
    }
  else
    outputlen2 = length2 - actg->cutTail - start2;

  if(isLargerThanTwin(ctg1))
    fprintf(fo, ">S%d_C%d_L%d_G%d", scafIndex, index_array[bal_ctg1], outputlen1, prevCtg->gapSeqLen);
  else
    fprintf(fo, ">S%d_C%d_L%d_G%d", scafIndex, index_array[ctg1], outputlen1, prevCtg->gapSeqLen);

  if(isLargerThanTwin(ctg2))
    fprintf(fo, "_C%d_L%d\n", index_array[bal_ctg2], outputlen2);
  else
    fprintf(fo, "_C%d_L%d\n", index_array[ctg2], outputlen2);

  if(contig_array[ctg1].seq)
    outputTightStr(fo, contig_array[ctg1].seq, start1, length1, outputlen1, 0, &column);
  else if(contig_array[bal_ctg1].seq)
    outputTightStr(fo, contig_array[bal_ctg1].seq, start1, length1, outputlen1, 1, &column);

  pt = (char *)darrayPut(gapSeqArray, prevCtg->gapSeqOffset);
  outputTightStrLowerCase(fo, pt, 0, prevCtg->gapSeqLen, prevCtg->gapSeqLen, 0, &column);

  if(contig_array[ctg2].seq)
    outputTightStr(fo, contig_array[ctg2].seq, start2, length2, outputlen2, 0, &column);
  else if(contig_array[bal_ctg2].seq)
    outputTightStr(fo, contig_array[bal_ctg2].seq, start2, length2, outputlen2, 1, &column);

  fprintf(fo, "\n");
}

static void outputGapSeq(FILE *fo, int index, STACK *ctgsStack, DARRAY *gapSeqArray)
{
  CTGinSCAF *actg, *prevCtg = NULL;
  stackRecover(ctgsStack);

  while((actg = stackPop(ctgsStack)) != NULL)
    {
      if(prevCtg && prevCtg->gapSeqLen > 0)
        output1gap(fo, index, prevCtg, actg, gapSeqArray);

      prevCtg = actg;

    }

}

static void outputScafSeq(FILE *fo, int index, STACK *ctgsStack, DARRAY *gapSeqArray)
{
  CTGinSCAF *actg, *prevCtg = NULL;
  unsigned int ctg, bal_ctg, length;
  int start, outputlen, gapN;
  char *pt;
  int column = 0;
  long long cvgSum = 0;
  int lenSum = 0;

  stackRecover(ctgsStack);

  while((actg = stackPop(ctgsStack)) != NULL)
    {
      if(!(contig_array[actg->ctgID].cvg > 0))
        continue;

      lenSum += contig_array[actg->ctgID].length;
      cvgSum += contig_array[actg->ctgID].length * contig_array[actg->ctgID].cvg;
    }

  if(lenSum > 0)
    fprintf(fo, ">scaffold%d %4.1f\n", index, (double)cvgSum / lenSum);
  else
    fprintf(fo, ">scaffold%d 0.0\n", index);

  stackRecover(ctgsStack);

  while((actg = stackPop(ctgsStack)) != NULL)
    {
      ctg = actg->ctgID;
      bal_ctg = getTwinCtg(ctg);
      length = contig_array[ctg].length + overlaplen;

      if(prevCtg && actg->scaftig_start)
        {
          gapN = actg->start - prevCtg->start - contig_array[prevCtg->ctgID].length;
          gapN = gapN > 0 ? gapN : 1;
          outputNs(fo, gapN, &column);
          //outputGapInfo(prevCtg->ctgID,ctg);
          Ncounter++;
        }

      if(!prevCtg)
        start = 0;
      else
        start = actg->cutHead;

      outputlen = length - start - actg->cutTail;

      if(contig_array[ctg].seq)
        outputTightStr(fo, contig_array[ctg].seq, start, length, outputlen, 0, &column);
      else if(contig_array[bal_ctg].seq)
        outputTightStr(fo, contig_array[bal_ctg].seq, start, length, outputlen, 1, &column);

      if(actg->gapSeqLen < 1)
        {
          prevCtg = actg;
          continue;
        }

      pt = (char *)darrayPut(gapSeqArray, actg->gapSeqOffset);
      outputTightStrLowerCase(fo, pt, 0, actg->gapSeqLen, actg->gapSeqLen, 0, &column);

      prevCtg = actg;
    }

  fprintf(fo, "\n");

}

static void fill1scaf(int index, STACK *ctgsStack, int thrdID);
static void check1scaf(int t, int thrdID)
{
  if(flagBuf[t])
    return;

  boolean late = 0;
  pthread_mutex_lock(&mutex);

  if(!flagBuf[t])
    {
      flagBuf[t] = 1;
      thrdNoBuf[t] = thrdID;
    }
  else
    late = 1;

  pthread_mutex_unlock(&mutex);

  if(late)
    return;

  counters[thrdID]++;
  fill1scaf(scafCounter + t + 1, ctgStackBuffer[t], thrdID);
}

static void fill1scaf(int index, STACK *ctgsStack, int thrdID)
{

  CTGinSCAF *actg, *prevCtg = NULL;
  READNEARBY *rdArray, *rdArray4gap, *rd;
  int numRd = 0, count, maxGLen = 0;
  unsigned int ctg, bal_ctg;
  STACK *rdStack;

  while((actg = stackPop(ctgsStack)) != NULL)
    {
      if(prevCtg)
        maxGLen = maxGLen < (actg->start - prevCtg->end) ? (actg->start - prevCtg->end) : maxGLen;

      ctg = actg->ctgID;
      bal_ctg = getTwinCtg(ctg);

      if(actg->mask)
        {
          prevCtg = actg;
          continue;
        }

      if(contig_array[ctg].closeReads)
        numRd += contig_array[ctg].closeReads->item_c;
      else if(contig_array[bal_ctg].closeReads)
        numRd += contig_array[bal_ctg].closeReads->item_c;

      prevCtg = actg;
    }

  if(numRd < 1)
    return;

  rdArray = (READNEARBY *)ckalloc(numRd * sizeof(READNEARBY));
  rdArray4gap = (READNEARBY *)ckalloc(numRd * sizeof(READNEARBY));
  //fprintf(stderr,"scaffold%d reads4gap %d\n",index,numRd);

  // collect reads appended to contigs in this scaffold
  int numRd2 = 0;
  stackRecover(ctgsStack);

  while((actg = stackPop(ctgsStack)) != NULL)
    {
      ctg = actg->ctgID;
      bal_ctg = getTwinCtg(ctg);

      if(actg->mask)
        continue;

      if(contig_array[ctg].closeReads)
        rdStack = contig_array[ctg].closeReads;
      else if(contig_array[bal_ctg].closeReads)
        rdStack = contig_array[bal_ctg].closeReads;
      else
        continue;

      stackBackup(rdStack);

      while((rd = (READNEARBY *)stackPop(rdStack)) != NULL)
        {
          rdArray[numRd2].len = rd->len;
          rdArray[numRd2].seqStarter = rd->seqStarter;

          if(isSmallerThanTwin(ctg))
            rdArray[numRd2++].dis = actg->start - overlaplen + rd->dis;
          else
            rdArray[numRd2++].dis = actg->start - overlaplen +
                                    contig_array[ctg].length - rd->len - rd->dis;
        }

      stackRecover(rdStack);
    }

  if(numRd2 != numRd)
    printf("##reads numbers doesn't match, %d vs %d when scaffold %d\n", numRd, numRd2, index);

  qsort(rdArray, numRd, sizeof(READNEARBY), cmp_reads);
  //fill gap one by one
  int gapStart, gapEnd;
  int numIn = 0;
  boolean flag;
  int buffer_size = maxReadLen > 100 ? maxReadLen : 100;
  int maxGSLen = maxGLen + GLDiff < 10 ? 10 : maxGLen + GLDiff;
  //fprintf(stderr,"maxGlen %d, maxGSlen %d\n",maxGLen,maxGSLen);

  char *seqGap = (char *)ckalloc(maxGSLen * sizeof(char)); // temp array for gap sequence
  Kmer *kmerCtg1 = (Kmer *)ckalloc(buffer_size * sizeof(Kmer));
  Kmer *kmerCtg2 = (Kmer *)ckalloc(buffer_size * sizeof(Kmer));
  char *seqCtg1 = (char *)ckalloc(buffer_size * sizeof(char));
  char *seqCtg2 = (char *)ckalloc(buffer_size * sizeof(char));
  prevCtg = NULL;
  stackRecover(ctgsStack);

  while((actg = stackPop(ctgsStack)) != NULL)
    {
      if(!prevCtg || !actg->scaftig_start)
        {
          prevCtg = actg;
          continue;
        }

      gapStart = prevCtg->end - 100;
      gapEnd = actg->start - overlaplen + 100;

      cutRdArray(rdArray, gapStart, gapEnd, &count, numRd, rdArray4gap);

      numIn += count;
      /*
      if(!count){
      	prevCtg = actg;
      	continue;
      }
      */
      int overlap;

      for(overlap = overlaplen; overlap > 14; overlap -= 2)
        {

          flag = localGraph(rdArray4gap, count, prevCtg, actg,
                            overlaplen, kmerCtg1, kmerCtg2, overlap, darrayBuf[thrdID],
                            seqCtg1, seqCtg2, seqGap);

          //free_kmerset(kmerSet);

          if(flag == 1)
            {
              /*
              fprintf(stderr,"Between ctg %d and %d, Found with %d\n",prevCtg->ctgID
              		,actg->ctgID,overlap);
              */
              break;
            }
        }

      /*
      if(count==0)
      	printf("Gap closed without reads\n");
      if(!flag)
      	fprintf(stderr,"Between ctg %d and %d, NO routes found\n",prevCtg->ctgID,actg->ctgID);
      */

      prevCtg = actg;
    }

  //fprintf(stderr,"____scaffold%d reads in gap %d\n",index,numIn);
  free((void *)seqGap);
  free((void *)kmerCtg1);
  free((void *)kmerCtg2);
  free((void *)seqCtg1);
  free((void *)seqCtg2);
  free((void *)rdArray);
  free((void *)rdArray4gap);
}

static void reverseStack(STACK *dStack, STACK *sStack)
{
  CTGinSCAF *actg, *ctgPt;
  emptyStack(dStack);

  while((actg = (CTGinSCAF *)stackPop(sStack)) != NULL)
    {
      ctgPt = (CTGinSCAF *)stackPush(dStack);
      ctgPt->ctgID = actg->ctgID;
      ctgPt->start = actg->start;
      ctgPt->end = actg->end;
      ctgPt->scaftig_start = actg->scaftig_start;
      ctgPt->mask = actg->mask;
      ctgPt->cutHead = actg->cutHead;
      ctgPt->cutTail = actg->cutTail;
      ctgPt->gapSeqLen = actg->gapSeqLen;
      ctgPt->gapSeqOffset = actg->gapSeqOffset;
    }

  stackBackup(dStack);
}

static Kmer tightStr2Kmer(char *tightStr, int start, int length, int revS)
{
  int i;
  Kmer word = 0;

  if(!revS)
    {
      if(start + overlaplen > length)
        {
          printf("tightStr2Kmer A: no enough bases for kmer\n");
          return word;
        }

      for(i = start; i < start + overlaplen; i++)
        {
          word <<= 2;
          word += getCharInTightString(tightStr, i);
        }
    }
  else
    {
      if(length - start - overlaplen < 0)
        {
          printf("tightStr2Kmer B: no enough bases for kmer\n");
          return word;
        }

      for(i = length - 1 - start; i > length - 1 - start - overlaplen; i--)
        {
          word <<= 2;
          word += int_comp(getCharInTightString(tightStr, i));
        }
    }

  return word;
}

static Kmer maxKmer()
{
  Kmer word = 0;
  int i;

  for(i = 0; i < overlaplen; i++)
    {
      word <<= 2;
      word += 0x3;
    }

  return word;
}
static int contigCatch(unsigned int prev_ctg, unsigned int ctg)
{
  if(contig_array[prev_ctg].length == 0 || contig_array[ctg].length == 0)
    return 0;

  Kmer kmerAtEnd, kmerAtStart;
  Kmer MaxKmer;
  unsigned int bal_ctg1 = getTwinCtg(prev_ctg);
  unsigned int bal_ctg2 = getTwinCtg(ctg);
  int i, start;
  int len1 = contig_array[prev_ctg].length + overlaplen;
  int len2 = contig_array[ctg].length + overlaplen;

  start = contig_array[prev_ctg].length;

  if(contig_array[prev_ctg].seq)
    kmerAtEnd = tightStr2Kmer(contig_array[prev_ctg].seq, start, len1, 0);
  else
    kmerAtEnd = tightStr2Kmer(contig_array[bal_ctg1].seq, start, len1, 1);

  start = 0;

  if(contig_array[ctg].seq)
    kmerAtStart = tightStr2Kmer(contig_array[ctg].seq, start, len2, 0);
  else
    kmerAtStart = tightStr2Kmer(contig_array[bal_ctg2].seq, start, len2, 1);

  MaxKmer = MAXKMER;

  for(i = 0; i < 10; i++)
    {
      if((kmerAtStart ^ kmerAtEnd) == 0)
        break;

      MaxKmer >>= 2;
      kmerAtEnd &= MaxKmer;
      kmerAtStart >>= 2;
    }

  if(i < 10)
    {
      return overlaplen - i;
    }
  else
    return 0;
}


static void initStackBuf(STACK **ctgStackBuffer, int scafBufSize)
{
  int i;

  for(i = 0; i < scafBufSize; i++)
    {
      flagBuf[i] = 1;
      ctgStackBuffer[i] = (STACK *)createStack(100, sizeof(CTGinSCAF));
    }
}
static void freeStackBuf(STACK **ctgStackBuffer, int scafBufSize)
{
  int i;

  for(i = 0; i < scafBufSize; i++)
    freeStack(ctgStackBuffer[i]);
}

static void threadRoutine(void *para)
{
  PARAMETER *prm;
  int i;

  prm = (PARAMETER *)para;

  //printf("%dth thread with threadID %d, hash_table %p\n",id,prm.threadID,prm.hash_table);
  while(1)
    {
      if(*(prm->selfSignal) == 1)
        {
          emptyDarray(darrayBuf[prm->threadID]);

          for(i = 0; i < scafInBuf; i++)
            check1scaf(i, prm->threadID);

          *(prm->selfSignal) = 0;
        }
      else if(*(prm->selfSignal) == 2)
        {
          *(prm->selfSignal) = 0;
          break;
        }

      usleep(1);
    }
}

static void creatThrds(pthread_t *threads, PARAMETER *paras)
{
  unsigned char i;
  int temp;

  for(i = 0; i < thrd_num; i++)
    {
      if((temp = pthread_create(&threads[i], NULL, (void *)threadRoutine, &(paras[i]))) != 0)
        {
          printf("create threads failed\n");
          exit(1);
        }
    }

  //printf("%d thread created\n...\n",thrd_num);
}

static void sendWorkSignal(unsigned char SIG, unsigned char *thrdSignals)
{
  int t;

  for(t = 0; t < thrd_num; t++)
    thrdSignals[t + 1] = SIG;

  while(1)
    {
      usleep(10);

      for(t = 0; t < thrd_num; t++)
        if(thrdSignals[t + 1])
          break;

      if(t == thrd_num)
        break;
    }
}

static void thread_wait(pthread_t *threads)
{
  int i;

  for(i = 0; i < thrd_num; i++)
    if(threads[i] != 0)
      pthread_join(threads[i], NULL);

}

static void outputSeqs(FILE *fo, FILE *fo2, int scafInBuf)
{
  int i, thrdID;

  for(i = 0; i < scafInBuf; i++)
    {
      thrdID = thrdNoBuf[i];
      outputScafSeq(fo, scafCounter + i + 1, ctgStackBuffer[i], darrayBuf[thrdID]);
      outputGapSeq(fo2, scafCounter + i + 1, ctgStackBuffer[i], darrayBuf[thrdID]);
    }
}

static void MaskContig(unsigned int ctg)
{
  contig_array[ctg].mask = 1;
  contig_array[getTwinCtg(ctg)].mask = 1;
}

static void MarkCtgOccu(unsigned int ctg)
{
  contig_array[ctg].flag = 1;
  contig_array[getTwinCtg(ctg)].flag = 1;
}

static void output_ctg(unsigned int ctg, FILE *fo)
{
  if(contig_array[ctg].length < 1)
    return;

  int len;
  unsigned int bal_ctg = getTwinCtg(ctg);

  len = contig_array[ctg].length + overlaplen;

  int col = 0;

  if(contig_array[ctg].seq)
    {
      fprintf(fo, ">C%d %4.1f\n", ctg, (double)contig_array[ctg].cvg);
      outputTightStr(fo, contig_array[ctg].seq, 0, len, len, 0, &col);
    }
  else if(contig_array[bal_ctg].seq)
    {
      fprintf(fo, ">C%d %4.1f\n", bal_ctg, (double)contig_array[ctg].cvg);
      outputTightStr(fo, contig_array[bal_ctg].seq, 0, len, len, 0, &col);
    }

  contig_array[ctg].flag = 1;
  contig_array[bal_ctg].flag = 1;
  fprintf(fo, "\n");
}

void prlReadsCloseGap(char *graphfile)
{
  //thrd_num=1;
  /*if(fillGap){
  	boolean flag;
  	//printf("\nStart to load reads for gap filling. %d length discrepancy is allowed\n",GLDiff);
  	//printf("...\n");
  	flag = loadReads4gap(graphfile);
  	if(!flag)
  		return;
  }*/

  if(orig2new)
    {
      convertIndex();
      orig2new = 0;
    }

  FILE *fp, *fo, *fo2;
  char line[1024];
  CTGinSCAF *actg;
  STACK *ctgStack, *aStack;
  int index = 0, offset = 0, counter, overallLen;
  int i, starter, prev_start, gapLen, catchable;
  unsigned int ctg, prev_ctg = 0;
  boolean IsPrevGap;
  pthread_t threads[thrd_num];
  unsigned char thrdSignal[thrd_num + 1];
  PARAMETER paras[thrd_num];

  for(ctg = 1; ctg <= num_ctg; ctg++)
    contig_array[ctg].flag = 0;

  MAXKMER = maxKmer();

  ctgStack = (STACK *)createStack(1000, sizeof(CTGinSCAF));

  sprintf(line, "%s.scaf_gap", graphfile);
  fp = ckopen(line, "r");
  sprintf(line, "%s.scafSeq", graphfile);
  fo = ckopen(line, "w");

  sprintf(line, "%s.gapSeq", graphfile);
  fo2 = ckopen(line, "w");

  pthread_mutex_init(&mutex, NULL);

  flagBuf = (boolean *)ckalloc(scafBufSize * sizeof(boolean));;
  thrdNoBuf = (unsigned char *)ckalloc(scafBufSize * sizeof(unsigned char));;
  memset(thrdNoBuf, 0, scafBufSize * sizeof(char));

  ctgStackBuffer = (STACK **)ckalloc(scafBufSize * sizeof(STACK *));
  initStackBuf(ctgStackBuffer, scafBufSize);

  darrayBuf = (DARRAY **)ckalloc(thrd_num * sizeof(DARRAY *));
  counters = (int *)ckalloc(thrd_num * sizeof(int));

  /*for(i=0;i<thrd_num;i++){
  	counters[i] = 0;
  	darrayBuf[i] = (DARRAY *)createDarray(100000,sizeof(char));
  	thrdSignal[i+1] = 0;
  	paras[i].threadID = i;
  	paras[i].mainSignal = &thrdSignal[0];
  	paras[i].selfSignal = &thrdSignal[i+1];
  }
  if(fillGap)
  	creatThrds(threads,paras);*/

  Ncounter = scafCounter = scafInBuf = allGaps = 0;

  while(fgets(line, sizeof(line), fp) != NULL)
    {
      if(line[0] == '>')
        {
          if(index)
            {
              aStack = ctgStackBuffer[scafInBuf];
              flagBuf[scafInBuf++] = 0;
              reverseStack(aStack, ctgStack);

              if(scafInBuf == scafBufSize)
                {
                  /*if(fillGap)
                  	sendWorkSignal(1,thrdSignal);*/

                  outputSeqs(fo, fo2, scafInBuf);
                  scafCounter += scafInBuf;
                  scafInBuf = 0;
                }

              //if(index%1000==0)
              //printf("Processed %d scaffolds\n",index);

            }

          //read next scaff
          emptyStack(ctgStack);
          IsPrevGap = offset = prev_ctg = 0;
          sscanf(line + 9, "%d %d %d", &index, &counter, &overallLen);
          continue;
        }

      if(line[0] == 'G') // gap appears
        {
          /*if(fillGap){
          	gapLen = procGap(line,ctgStack);
          	IsPrevGap = 1;
          }*/
          continue;
        }

      if(line[0] >= '0' && line[0] <= '9') // a contig line
        {
          sscanf(line, "%d %d", &ctg, &starter);
          actg = (CTGinSCAF *)stackPush(ctgStack);
          actg->ctgID = ctg;

          if(contig_array[ctg].flag)
            MaskContig(ctg);
          else
            MarkCtgOccu(ctg);

          initiateCtgInScaf(actg);

          if(!prev_ctg)
            actg->cutHead = 0;
          else if(!IsPrevGap)
            allGaps++;

          if(!IsPrevGap)
            {
              if(prev_ctg && (starter - prev_start - (int)contig_array[prev_ctg].length)
                  < ((int)overlaplen * 4))
                {
                  /*
                  if(fillGap)
                  	catchable = contigCatch(prev_ctg,ctg);
                  else
                  */
                  catchable = 0;

                  if(catchable)   // prev_ctg and ctg overlap **bp
                    {
                      allGaps--;
                      /*
                      if(isLargerThanTwin(prev_ctg))
                      	fprintf(stderr,"%d #######  by_overlap\n",getTwinCtg(prev_ctg));
                      else
                      	fprintf(stderr,"%d ####### by_overlap\n",prev_ctg);
                      */
                      actg->scaftig_start = 0;
                      actg->cutHead = catchable;
                      offset += - (starter - prev_start - contig_array[prev_ctg].length) +
                                (overlaplen - catchable);
                    }
                  else
                    actg->scaftig_start = 1;

                }
              else
                actg->scaftig_start = 1;
            }
          else
            {
              offset += - (starter - prev_start - contig_array[prev_ctg].length) + gapLen;
              actg->scaftig_start = 0;
            }

          actg->start = starter + offset;
          actg->end = actg->start + contig_array[ctg].length - 1;
          actg->mask = contig_array[ctg].mask;
          IsPrevGap = 0;
          prev_ctg = ctg;
          prev_start = starter;
        }
    }

  if(index)
    {
      aStack = ctgStackBuffer[scafInBuf];
      flagBuf[scafInBuf++] = 0;
      reverseStack(aStack, ctgStack);

      if(fillGap)
        sendWorkSignal(1, thrdSignal);

      outputSeqs(fo, fo2, scafInBuf);
    }

  /*if(fillGap){
  	sendWorkSignal(2,thrdSignal);
  	thread_wait(threads);
  }*/
  for(ctg = 1; ctg <= num_ctg; ctg++)
    {
      if((contig_array[ctg].length + overlaplen) < 100 ||
          contig_array[ctg].flag)
        continue;

      output_ctg(ctg, fo);

    }

  //printf("Done with %d scaffolds, %d gaps finished, %d gaps overall\n",index,allGaps-Ncounter,allGaps);
  //printf("scaffolds outputted : %d.\n",index);
  index = 0;

  for(i = 0; i < thrd_num; i++)
    {
      freeDarray(darrayBuf[i]);
      index += counters[i];
    }

  if(fillGap)
    //printf("Threads processed %d scaffolds\n",index);
    free((void *)darrayBuf);

  if(readSeqInGap)
    freeDarray(readSeqInGap);

  fclose(fp);
  fclose(fo);
  fclose(fo2);
  freeStack(ctgStack);
  freeStackBuf(ctgStackBuffer, scafBufSize);
  free((void *)flagBuf);
  free((void *)thrdNoBuf);
  free((void *)ctgStackBuffer);
}