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

#define CTGendLen 35  // shouldn't larger than max_read_len
#define UPlimit 5000
#define MaxRouteNum 10

static Kmer pubKmer = 0x1b4d65165b;

static void kmerSet_mark(KmerSet *set);
static void trace4Repeat(Kmer currW, int steps, int min, int max, int *num_route,
                         KmerSet *kset, Kmer kmerDest, int overlap, Kmer WORDF,
                         int *traceCounter, int maxRoute, kmer_t **soFarNode, short *multiOccu1, short *multiOccu2,
                         int *routeLens, char **foundRoutes, char *soFarSeq,
                         long long *soFarLinks, double *avgLinks);

static Kmer prevKmerLocal(Kmer next, char ch, int overlap)
{
  Kmer word = next;
  word >>= 2;
  word += ((Kmer)ch) << 2 * (overlap - 1);
  return word;
}
static Kmer nextKmerLocal(Kmer prev, char ch, Kmer WordFilter)
{
  Kmer word = prev;
  word <<= 2;
  word &= WordFilter;
  word += ch;
  return word;
}
static void singleKmer(int t, KmerSet *kset, int flag, Kmer *kmerBuffer, char *prevcBuffer, char *nextcBuffer)
{
  kmer_t *pos;

  put_kmerset(kset, kmerBuffer[t], prevcBuffer[t], nextcBuffer[t], &pos);

  if(pos->inEdge == flag)
    return;
  else if(pos->inEdge == 0)
    pos->inEdge = flag;
  else if(pos->inEdge == 1 && flag == 2)
    pos->inEdge = 3;
  else if(pos->inEdge == 2 && flag == 1)
    pos->inEdge = 3;

}

static void putKmer2DBgraph(KmerSet *kset, int flag, int kmer_c, Kmer *kmerBuffer, char *prevcBuffer, char *nextcBuffer)
{
  int t;

  for(t = 0; t < kmer_c; t++)
    singleKmer(t, kset, flag, kmerBuffer, prevcBuffer, nextcBuffer);

}

static void getSeqFromRead(READNEARBY read, char *src_seq)
{
  int len_seq = read.len;
  int j;
  char *tightSeq = (char *)darrayGet(readSeqInGap, read.seqStarter);

  for(j = 0; j < len_seq; j++)
    src_seq[j] = getCharInTightString(tightSeq, j);
}

static void chopKmer4Ctg(Kmer *kmerCtg, int lenCtg, int overlap, char *src_seq, Kmer WORDF)
{
  int index, j;
  Kmer word = 0;

  for (index = 0; index < overlap; index++)
    {
      word <<= 2;
      word += src_seq[index];
    }

  index = 0;
  kmerCtg[index++] = word;

  for(j = 1; j <= lenCtg - overlap; j ++)
    {
      word = nextKmerLocal(word, src_seq[j - 1 + overlap], WORDF);
      kmerCtg[index++] = word;
    }
}

static void chopKmer4read(int len_seq, int overlap, char *src_seq, char *bal_seq,
                          Kmer *kmerBuffer, char *prevcBuffer, char *nextcBuffer, int *kmer_c, Kmer WORDF)
{
  int j, bal_j;
  Kmer word, bal_word;
  int index;
  char InvalidCh = 4;

  if(len_seq < overlap + 1)
    {
      *kmer_c = 0;
      return;
    }

  word = 0;

  for (index = 0; index < overlap; index++)
    {
      word <<= 2;
      word += src_seq[index];
    }

  reverseComplementSeq(src_seq, len_seq, bal_seq);

  // complementary node
  bal_word = reverseComplement(word, overlap);
  bal_j = len_seq - 0 - overlap; //  0;
  index = 0;

  if(word < bal_word)
    {
      kmerBuffer[index] = word;
      prevcBuffer[index] = InvalidCh;
      nextcBuffer[index++] = src_seq[0 + overlap];
    }
  else
    {
      kmerBuffer[index] = bal_word;
      prevcBuffer[index] = bal_seq[bal_j - 1];
      nextcBuffer[index++] = InvalidCh;
    }

  for(j = 1; j <= len_seq - overlap; j ++)
    {
      word = nextKmerLocal(word, src_seq[j - 1 + overlap], WORDF);
      bal_j = len_seq - j - overlap; //  j;
      bal_word = prevKmerLocal(bal_word, bal_seq[bal_j], overlap);

      if(word < bal_word)
        {
          kmerBuffer[index] = word;
          prevcBuffer[index] = src_seq[j - 1];

          if(j < len_seq - overlap)
            nextcBuffer[index++] = src_seq[j + overlap];
          else
            nextcBuffer[index++] = InvalidCh;

          //printf("%dth: %p with %p\n",kmer_c-1,word,hashBanBuffer[kmer_c-1]);
        }
      else
        {
          // complementary node
          kmerBuffer[index] = bal_word;

          if(bal_j > 0)
            prevcBuffer[index] = bal_seq[bal_j - 1];
          else
            prevcBuffer[index] = InvalidCh;

          nextcBuffer[index++] = bal_seq[bal_j + overlap];
          //printf("%dth: %p with %p\n",kmer_c-1,bal_word,hashBanBuffer[kmer_c-1]);
        }
    }

  *kmer_c = index;
}

static void headTightStr(char *tightStr, int length, int start, int headLen, int revS, char *src_seq)
{
  int i, index = 0;

  if(!revS)
    {
      for(i = start; i < start + headLen; i++)
        src_seq[index++] = getCharInTightString(tightStr, i);
    }
  else
    {
      for(i = length - 1 - start; i >= length - headLen - start; i--)
        src_seq[index++] = int_comp(getCharInTightString(tightStr, i));
    }
}

static int getSeqFromCtg(CTGinSCAF *ctg, boolean fromHead, unsigned int len, int originOverlap, char *src_seq)
{
  unsigned int ctgId = ctg->ctgID;
  unsigned int bal_ctg = getTwinCtg(ctgId);

  if(contig_array[ctgId].length < 1)
    return 0;

  unsigned int length = contig_array[ctgId].length + originOverlap;

  len = len < length ? len : length;

  if(fromHead)
    {
      if(contig_array[ctgId].seq)
        headTightStr(contig_array[ctgId].seq, length, 0, len, 0, src_seq);
      else
        headTightStr(contig_array[bal_ctg].seq, length, 0, len, 1, src_seq);
    }
  else
    {
      if(contig_array[ctgId].seq)
        headTightStr(contig_array[ctgId].seq, length, length - len, len, 0, src_seq);
      else
        headTightStr(contig_array[bal_ctg].seq, length, length - len, len, 1, src_seq);
    }

  return len;
}


static KmerSet *readsInGap2DBgraph(READNEARBY *rdArray, int num, CTGinSCAF *ctg1, CTGinSCAF *ctg2, int originOverlap,
                                   Kmer *kmerCtg1, Kmer *kmerCtg2, int overlap, Kmer WordFilter)
{
  int kmer_c;
  Kmer *kmerBuffer;
  char *nextcBuffer, *prevcBuffer;
  int i;
  int buffer_size = maxReadLen > CTGendLen ? maxReadLen : CTGendLen;
  KmerSet *kmerS = NULL;
  int lenCtg1;
  int lenCtg2;
  char *bal_seq;
  char *src_seq;

  src_seq = (char *)ckalloc(buffer_size * sizeof(char));
  bal_seq = (char *)ckalloc(buffer_size * sizeof(char));

  kmerBuffer = (Kmer *)ckalloc(buffer_size * sizeof(Kmer));
  prevcBuffer = (char *)ckalloc(buffer_size * sizeof(char));
  nextcBuffer = (char *)ckalloc(buffer_size * sizeof(char));

  kmerS = init_kmerset(1024, 0.77f);

  for(i = 0; i < num; i++)
    {
      getSeqFromRead(rdArray[i], src_seq);
      chopKmer4read(rdArray[i].len, overlap, src_seq, bal_seq,
                    kmerBuffer, prevcBuffer, nextcBuffer, &kmer_c, WordFilter);
      putKmer2DBgraph(kmerS, 0, kmer_c, kmerBuffer, prevcBuffer, nextcBuffer);
    }

  lenCtg1 = getSeqFromCtg(ctg1, 0, CTGendLen, originOverlap, src_seq);
  chopKmer4Ctg(kmerCtg1, lenCtg1, overlap, src_seq, WordFilter);
  chopKmer4read(lenCtg1, overlap, src_seq, bal_seq,
                kmerBuffer, prevcBuffer, nextcBuffer, &kmer_c, WordFilter);
  putKmer2DBgraph(kmerS, 1, kmer_c, kmerBuffer, prevcBuffer, nextcBuffer);

  lenCtg2 = getSeqFromCtg(ctg2, 1, CTGendLen, originOverlap, src_seq);
  chopKmer4Ctg(kmerCtg2, lenCtg2, overlap, src_seq, WordFilter);
  chopKmer4read(lenCtg2, overlap, src_seq, bal_seq,
                kmerBuffer, prevcBuffer, nextcBuffer, &kmer_c, WordFilter);
  putKmer2DBgraph(kmerS, 2, kmer_c, kmerBuffer, prevcBuffer, nextcBuffer);
  /*
  	if(ctg1->ctgID==3733&&ctg2->ctgID==3067){
  		for(i=0;i<lenCtg2;i++)
  			printf("%c",int2base((int)src_seq[i]));
  		printf("\n");
  	}
  */
  //printf("sequence length chop from contigs on both sides: %d %d\n",lenCtg1,lenCtg2);
  //kmerSet_deLoop(kmerS,WordFilter);
  kmerSet_mark(kmerS);
  free((void *)src_seq);
  free((void *)bal_seq);
  free((void *)kmerBuffer);
  free((void *)nextcBuffer);
  free((void *)prevcBuffer);

  fflush(stdout);

  return kmerS;
}

static void printKmer(FILE *fo, Kmer kmer, int overlap)
{
  int i;
  char kmerSeq[32], ch;

  for(i = overlap - 1; i >= 0; i--)
    {
      ch = kmer & 3;
      kmer >>= 2;
      kmerSeq[i] = ch;
    }

  for(i = 0; i < overlap; i++)
    fprintf(fo, "%c", int2base((int)kmerSeq[i]));
}

static void kmerSet_mark(KmerSet *set)
{
  int i, in_num, out_num, cvgSingle;
  kmer_t *rs;
  long long counter = 0, linear = 0;
  Kmer word;

  set->iter_ptr = 0;

  while(set->iter_ptr < set->size)
    {
      if(!is_kmer_entity_null(set->flags, set->iter_ptr))
        {
          in_num = out_num = 0;
          rs = set->array + set->iter_ptr;
          word = rs->seq;

          for(i = 0; i < 4; i++)
            {
              cvgSingle = get_kmer_left_cov(*rs, i);

              if(cvgSingle > 0)
                {
                  in_num++;
                }

              cvgSingle = get_kmer_right_cov(*rs, i);

              if(cvgSingle > 0)
                {
                  out_num++;
                }
            }

          if(rs->single)
            {
              counter++;
            }

          if(in_num == 1 && out_num == 1)
            {
              rs->linear = 1;
              linear++;
            }
        }

      set->iter_ptr ++;
    }

  //printf("Allocated %ld node, %ld single nodes, %ld linear\n",(long)count_kmerset(set),counter,linear);
}

static kmer_t *searchNode(Kmer word, KmerSet *kset, int overlap)
{
  Kmer bal_word = reverseComplement(word, overlap);
  kmer_t *node;
  boolean found;

  if(word < bal_word)
    found = search_kmerset(kset, word, &node);
  else
    found = search_kmerset(kset, bal_word, &node);

  if(found)
    return node;
  else
    return NULL;
}

static int searchKmerOnCtg(Kmer currW, Kmer *kmerDest, int num)
{
  int i;

  for(i = 0; i < num; i++)
    {
      if(currW == kmerDest[i])
        {
          return i;
        }
    }

  return -1;
}

// pick on from n items randomly
static int nPick1(int *array, int n)
{
  int m, i;
  m = n - 1; //(int)(drand48()*n);
  int value = array[m];

  for(i = m; i < n - 1; i++)
    array[i] = array[i + 1];

  return value;
}

static void traceAlongDBgraph(Kmer currW, int steps, int min, int max, int *num_route,
                              KmerSet *kset, Kmer *kmerDest, int num, int overlap, Kmer WORDF,
                              char **foundRoutes, int *routeEndOnCtg2, int *routeLens, char *soFarSeq,
                              int *traceCounter, int maxRoute, kmer_t **soFarNode, boolean *multiOccu,
                              long long *soFarLinks, double *avgLinks)
{
  (*traceCounter)++;

  if(*traceCounter > UPlimit)
    {
      /*
      if(overlap==19&&kmerDest[0]==pubKmer)
      	printf("UPlimit\n");
      */
      return;
    }

  if(steps > max || *num_route >= maxRoute)
    {
      /*
      if(overlap==19&&kmerDest[0]==pubKmer)
      	printf("max steps/maxRoute\n");
      */
      return;
    }

  Kmer word = reverseComplement(currW, overlap);
  boolean isSmaller = currW < word;
  int i;
  char ch;
  unsigned char links;

  if(isSmaller)
    word = currW;

  kmer_t *node;
  boolean found = search_kmerset(kset, word, &node);

  if(!found)
    {
      printf("Trace: can't find kmer %llx (rc %llx, input %llx) at step %d\n", word,
             reverseComplement(word, overlap), currW, steps);
      return;
    }

  if(node->twin > 1)
    return;

  if(soFarNode)
    soFarNode[steps] = node;

  if(steps > 0)
    soFarSeq[steps - 1] = currW & 0x03;

  int index, end;
  int linkCounter = *soFarLinks;

  if(steps >= min && node->inEdge > 1 && (end = searchKmerOnCtg(currW, kmerDest, num)) >= 0)
    {
      index = *num_route;

      if(steps > 0)
        avgLinks[index] = (double)linkCounter / steps;
      else
        avgLinks[index] = 0;

      //find node that appears more than once in the path
      multiOccu[index] = 0;

      for(i = 0; i < steps + 1; i++)
        soFarNode[i]->deleted = 0;

      for(i = 0; i < steps + 1; i++)
        {
          if(soFarNode[i]->deleted)
            {
              multiOccu[index] = 1;
              break;
            }

          soFarNode[i]->deleted = 1;
        }

      routeEndOnCtg2[index] = end;
      routeLens[index] = steps;
      char *array = foundRoutes[index];

      for(i = 0; i < steps; i++)
        array[i] = soFarSeq[i];

      if(i < max)
        array[i] = 4;  //indicate the end of the sequence

      *num_route = ++index;
      return;
    }

  steps++;

  if(isSmaller)
    {
      int array[] = {0, 1, 2, 3};

      for(i = 4; i > 0; i--)
        {
          ch = nPick1(array, i);
          links = get_kmer_right_cov(*node, ch);

          if(!links)
            continue;

          *soFarLinks = linkCounter + links;
          word = nextKmerLocal(currW, ch, WORDF);
          traceAlongDBgraph(word, steps, min, max, num_route,
                            kset, kmerDest, num, overlap, WORDF,
                            foundRoutes, routeEndOnCtg2, routeLens, soFarSeq,
                            traceCounter, maxRoute, soFarNode, multiOccu,
                            soFarLinks, avgLinks);
        }
    }
  else
    {
      int array[] = {0, 1, 2, 3};

      for(i = 4; i > 0; i--)
        {
          ch = nPick1(array, i);
          links = get_kmer_left_cov(*node, ch);

          if(!links)
            continue;

          *soFarLinks = linkCounter + links;
          word = nextKmerLocal(currW, int_comp(ch), WORDF);
          traceAlongDBgraph(word, steps, min, max, num_route,
                            kset, kmerDest, num, overlap, WORDF,
                            foundRoutes, routeEndOnCtg2, routeLens, soFarSeq,
                            traceCounter, maxRoute, soFarNode, multiOccu,
                            soFarLinks, avgLinks);
        }
    }
}

static int searchFgap(KmerSet *kset, CTGinSCAF *ctg1, CTGinSCAF *ctg2, Kmer *kmerCtg1,
                      Kmer *kmerCtg2, unsigned int origOverlap, int overlap, DARRAY *gapSeqArray,
                      int len1, int len2, Kmer WordFilter, int *offset1, int *offset2, char *seqGap, int *cut1, int *cut2)
{

  int i;
  int ret = 0;
  kmer_t *node, **soFarNode;
  int num_route;
  int gapLen = ctg2->start - ctg1->end - origOverlap + overlap;
  int min = gapLen - GLDiff > 0 ? gapLen - GLDiff : 0; //0531
  int max = gapLen + GLDiff < 10 ? 10 : gapLen + GLDiff;
  char **foundRoutes;
  char *soFarSeq;
  int traceCounter;
  int *routeEndOnCtg2;
  int *routeLens;
  boolean *multiOccu;
  long long soFarLinks;
  double *avgLinks;

  //mask linear internal linear kmer on contig1 end
  routeEndOnCtg2 = (int *)ckalloc(MaxRouteNum * sizeof(int));
  routeLens = (int *)ckalloc(MaxRouteNum * sizeof(int));
  multiOccu = (boolean *)ckalloc(MaxRouteNum * sizeof(boolean));
  short *MULTI1 = (short *)ckalloc(MaxRouteNum * sizeof(short));
  short *MULTI2 = (short *)ckalloc(MaxRouteNum * sizeof(short));
  soFarSeq = (char *)ckalloc(max * sizeof(char));
  soFarNode = (kmer_t **)ckalloc((max + 1) * sizeof(kmer_t *));
  foundRoutes = (char **)ckalloc(MaxRouteNum * sizeof(char *));;
  avgLinks = (double *)ckalloc(MaxRouteNum * sizeof(double));;

  for(i = 0; i < MaxRouteNum; i++)
    foundRoutes[i] = (char *)ckalloc(max * sizeof(char));

  for(i = len1 - 1; i >= 0; i--)
    {

      num_route = traceCounter = soFarLinks = 0;
      int steps = 0;
      traceAlongDBgraph(kmerCtg1[i], steps, min, max, &num_route,
                        kset, kmerCtg2, len2, overlap, WordFilter,
                        foundRoutes, routeEndOnCtg2, routeLens, soFarSeq,
                        &traceCounter, MaxRouteNum, soFarNode, multiOccu,
                        &soFarLinks, avgLinks);

      if(num_route > 0)
        {
          int m, minEnd = routeEndOnCtg2[0];

          for(m = 0; m < num_route; m++)
            {
              if(routeLens[m] < 0)
                continue;

              if(routeEndOnCtg2[m] < minEnd)
                minEnd = routeEndOnCtg2[m];
            }

          /* else if(minFreq>1){
          for(m=0;m<num_route;m++){
          	if(routeEndOnCtg2[m]!=minEnd)
          		continue;
          	for(j=0;j<max;j++){
          		if(foundRoutes[m][j]>3)
          			break;
          		printf("%c",int2base((int)foundRoutes[m][j]));
          	}
          	printf(": %4.2f\n",avgLinks[m]);
          }
          }   */

          num_route = traceCounter = soFarLinks = 0;
          steps = 0;
          trace4Repeat(kmerCtg1[i], steps, min, max, &num_route,
                       kset, kmerCtg2[minEnd], overlap, WordFilter,
                       &traceCounter, MaxRouteNum, soFarNode, MULTI1, MULTI2,
                       routeLens, foundRoutes, soFarSeq, &soFarLinks, avgLinks);
          int j, best = 0;
          int maxLen = routeLens[0];
          double maxLink = avgLinks[0];
          char *pt;
          boolean repeat = 0, sameLen = 1;
          int leftMost = max, rightMost = max;

          if(num_route < 1)
            {
              fprintf(stderr, "After trace4Repeat: non route was found\n");
              continue;
            }

          if(num_route > 1)
            {
              // if multi paths are found, we check on the repeatative occurrences and links/length
              for(m = 0; m < num_route; m++)
                {
                  if(routeLens[m] < 0)
                    continue;

                  if(MULTI1[m] >= 0 && MULTI2[m] >= 0)
                    {
                      repeat = 1;
                      leftMost = leftMost > MULTI1[m] ? MULTI1[m] : leftMost;
                      rightMost = rightMost > MULTI2[m] ? MULTI2[m] : rightMost;
                    }

                  if(routeLens[m] != maxLen)
                    sameLen = 0;

                  if(routeLens[m] < maxLen)
                    maxLen = routeLens[m];

                  if(avgLinks[m] > maxLink)
                    {
                      maxLink = avgLinks[m];
                      best = m;
                    }
                }
            }

          if(repeat)
            {
              *offset1 = *offset2 = *cut1 = *cut2 = 0;
              int index = 0;
              char ch;

              for(j = 0; j < leftMost; j++)
                {
                  if(routeLens[0] < j + overlap + 1)
                    break;
                  else
                    ch = foundRoutes[0][j];

                  for(m = 1; m < num_route; m++)
                    {
                      if(routeLens[m] < 0)
                        continue;

                      if(ch != foundRoutes[m][j])
                        break;
                    }

                  if(m == num_route)
                    seqGap[index++] = ch;
                  else break;
                }

              *offset1 = index;
              index = 0;

              for(j = 0; j < rightMost; j++)
                {
                  if(routeLens[0] - overlap - 1 < j)
                    break;
                  else
                    ch = foundRoutes[0][routeLens[0] - overlap - 1 - j];

                  for(m = 1; m < num_route; m++)
                    {
                      if(routeLens[m] < 0)
                        continue;

                      if(ch != foundRoutes[m][routeLens[m] - overlap - 1 - j])
                        break;
                    }

                  if(m == num_route)
                    index++;
                  else break;
                }

              *offset2 = index;

              for(j = 0; j < *offset2; j++)
                seqGap[*offset1 + *offset2 - 1 - j] = foundRoutes[0][routeLens[0] - overlap - 1 - j];

              if(*offset1 > 0 || *offset2 > 0)
                {
                  *cut1 = len1 - i - 1;
                  *cut2 = minEnd;

                  //fprintf(stderr,"\n");
                  for(m = 0; m < num_route; m++)
                    {
                      for(j = 0; j < max; j++)
                        {
                          if(foundRoutes[m][j] > 3)
                            break;

                          //fprintf(stderr,"%c",int2base((int)foundRoutes[m][j]));
                        }

                      //fprintf(stderr,": %4.2f\n",avgLinks[m]);
                    }

                  /*
                  fprintf(stderr,">Gap (%d + %d) (%d + %d)\n",*offset1,*offset2,*cut1,*cut2);
                  for(index=0;index<*offset1+*offset2;index++)
                  	fprintf(stderr,"%c",int2base(seqGap[index]));
                  fprintf(stderr,"\n");  */
                }

              ret = 3;
              break;
            }

          if(overlap + (len1 - i - 1) + minEnd - routeLens[best] > (int)origOverlap)
            continue;

          ctg1->gapSeqOffset = gapSeqArray->item_c;
          ctg1->gapSeqLen = routeLens[best];

          if(!darrayPut(gapSeqArray, ctg1->gapSeqOffset + maxLen / 4))
            continue;

          pt = (char *)darrayPut(gapSeqArray, ctg1->gapSeqOffset);

          /*
          printKmer(stderr,kmerCtg1[i],overlap);
          fprintf(stderr,"-");
          */
          for(j = 0; j < max; j++)
            {
              if(foundRoutes[best][j] > 3)
                break;

              writeChar2tightString(foundRoutes[best][j], pt, j);
              //fprintf(stderr,"%c",int2base((int)foundRoutes[best][j]));
            }

          //fprintf(stderr,": GAPSEQ %d + %d, avglink %4.2f\n",len1-i-1,minEnd,avgLinks[best]);
          ctg1->cutTail = len1 - i - 1;
          ctg2->cutHead = overlap + minEnd;
          ctg2->scaftig_start = 0;

          ret = 1;
          break;
          /* }if(num_route>1){
          	ret = 2;
          	break;    */
        }
      else    //mark node which leads to dead end
        {
          node = searchNode(kmerCtg1[i], kset, overlap);

          if(node)
            node->twin = 2;
        }

    }

  for(i = 0; i < MaxRouteNum; i++)
    free((void *)foundRoutes[i]);

  free((void *)soFarSeq);
  free((void *)soFarNode);
  free((void *)multiOccu);
  free((void *)MULTI1);
  free((void *)MULTI2);
  free((void *)foundRoutes);
  free((void *)routeEndOnCtg2);
  free((void *)routeLens);

  return ret;
}

static void trace4Repeat(Kmer currW, int steps, int min, int max, int *num_route,
                         KmerSet *kset, Kmer kmerDest, int overlap, Kmer WORDF,
                         int *traceCounter, int maxRoute, kmer_t **soFarNode, short *multiOccu1, short *multiOccu2,
                         int *routeLens, char **foundRoutes, char *soFarSeq,
                         long long *soFarLinks, double *avgLinks)
{
  (*traceCounter)++;

  if(*traceCounter > UPlimit)
    return;

  if(steps > max || *num_route >= maxRoute)
    return;

  Kmer word = reverseComplement(currW, overlap);
  boolean isSmaller = currW < word;
  char ch;
  unsigned char links;
  int index, i;

  if(isSmaller)
    word = currW;

  kmer_t *node;
  boolean found = search_kmerset(kset, word, &node);

  if(!found)
    {
      printf("Trace: can't find kmer %llx (rc %llx, input %llx) at step %d\n", word,
             reverseComplement(word, overlap), currW, steps);
      return;
    }

  if(soFarNode)
    soFarNode[steps] = node;

  if(soFarSeq && steps > 0)
    soFarSeq[steps - 1] = currW & 0x03;

  int linkCounter;

  if(soFarLinks)
    linkCounter = *soFarLinks;

  if(steps >= min && currW == kmerDest)
    {
      index = *num_route;

      if(avgLinks && steps > 0)
        avgLinks[index] = (double)linkCounter / steps;
      else if(avgLinks)
        avgLinks[index] = 0;

      //find node that appears more than once in the path
      if(multiOccu1 && multiOccu2)
        {
          for(i = 0; i < steps + 1; i++)
            soFarNode[i]->deleted = 0;

          int rightMost = 0;
          boolean MULTI = 0;

          for(i = 0; i < steps + 1; i++)
            {
              if(soFarNode[i]->deleted)
                {
                  rightMost = rightMost < i - 1 ? i - 1 : rightMost;
                  MULTI = 1;
                }

              soFarNode[i]->deleted = 1;
            }

          if(!MULTI)
            multiOccu1[index] = multiOccu2[index] = -1;
          else
            {
              multiOccu2[index] = steps - 2 - rightMost < 0 ? 0 : steps - 2 - rightMost; //[0 steps-2]

              for(i = 0; i < steps + 1; i++)
                soFarNode[i]->deleted = 0;

              int leftMost = steps - 2;

              for(i = steps; i >= 0; i--)
                {
                  if(soFarNode[i]->deleted)
                    leftMost = leftMost > i - 1 ? i - 1 : leftMost;

                  soFarNode[i]->deleted = 1;
                }

              multiOccu1[index] = leftMost < 0 ? 0 : leftMost; //[0 steps-2]
            }
        }

      if(routeLens)
        routeLens[index] = steps;

      if(soFarSeq)
        {
          char *array = foundRoutes[index];

          for(i = 0; i < steps; i++)
            array[i] = soFarSeq[i];

          if(i < max)
            array[i] = 4;  //indicate the end of the sequence
        }

      *num_route = ++index;
    }

  steps++;

  if(isSmaller)
    {
      int array[] = {0, 1, 2, 3};

      for(i = 4; i > 0; i--)
        {
          ch = nPick1(array, i);
          links = get_kmer_right_cov(*node, ch);

          if(!links)
            continue;

          if(soFarLinks)
            *soFarLinks = linkCounter + links;

          word = nextKmerLocal(currW, ch, WORDF);
          trace4Repeat(word, steps, min, max, num_route,
                       kset, kmerDest, overlap, WORDF, traceCounter, maxRoute, soFarNode,
                       multiOccu1, multiOccu2, routeLens, foundRoutes, soFarSeq,
                       soFarLinks, avgLinks);
        }
    }
  else
    {
      int array[] = {0, 1, 2, 3};

      for(i = 4; i > 0; i--)
        {
          ch = nPick1(array, i);
          links = get_kmer_left_cov(*node, ch);

          if(!links)
            continue;

          if(soFarLinks)
            *soFarLinks = linkCounter + links;

          word = nextKmerLocal(currW, int_comp(ch), WORDF);
          trace4Repeat(word, steps, min, max, num_route,
                       kset, kmerDest, overlap, WORDF, traceCounter, maxRoute, soFarNode,
                       multiOccu1, multiOccu2, routeLens, foundRoutes, soFarSeq,
                       soFarLinks, avgLinks);
        }
    }
}

//found repeat node on contig ends
static void maskRepeatNode(KmerSet *kset, Kmer *kmerCtg1,
                           Kmer *kmerCtg2, int overlap,
                           int len1, int len2, int max, Kmer WordFilter)
{
  int i;
  int num_route, steps;
  int min = 1, maxRoute = 1;
  int traceCounter;
  Kmer word, bal_word;
  kmer_t *node;
  boolean found;
  int counter = 0;

  for(i = 0; i < len1; i++)
    {
      word = kmerCtg1[i];
      bal_word = reverseComplement(word, overlap);

      if(word > bal_word)
        word = bal_word;

      found = search_kmerset(kset, word, &node);

      if(!found || node->linear)
        {
          //printf("Found no node for kmer %llx\n",word);
          continue;
        }

      num_route = traceCounter = 0;
      steps = 0;
      trace4Repeat(word, steps, min, max, &num_route,
                   kset, word, overlap, WordFilter,
                   &traceCounter, maxRoute, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);

      if(num_route < 1)
        continue;

      counter++;
      node->checked = 1;
    }

  for(i = 0; i < len2; i++)
    {
      word = kmerCtg2[i];
      bal_word = reverseComplement(word, overlap);

      if(word > bal_word)
        word = bal_word;

      found = search_kmerset(kset, word, &node);

      if(!found || node->linear)
        {
          //printf("Found no node for kmer %llx\n",word);
          continue;
        }

      num_route = traceCounter = 0;
      steps = 0;
      trace4Repeat(word, steps, min, max, &num_route,
                   kset, word, overlap, WordFilter,
                   &traceCounter, maxRoute, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);

      if(num_route < 1)
        continue;

      counter++;
      node->checked = 1;
    }

  //printf("MR: %d(%d)\n",counter,len1+len2);
}

/*
static boolean chopReadFillGap(int len_seq,int overlap,char *src_seq, char *bal_seq,
					KmerSet *kset,Kmer WORDF,int *start,int *end,boolean *bal,
					Kmer *KmerCtg1,int len1,Kmer *KmerCtg2,int len2,int *index1,int *index2)
{
	int index,j=0,bal_j;
	Kmer word,bal_word;
	int flag=0,bal_flag=0;
	int ctg1start,bal_ctg1start,ctg2end,bal_ctg2end;
	int seqStart,bal_start,seqEnd,bal_end;
	kmer_t *node;
	boolean found;

	if(len_seq<overlap+1){
		return 0;
	}
	word = 0;
	for (index = 0;index<overlap;index++){
		word <<= 2;
		word += src_seq[index];
	}
	reverseComplementSeq(src_seq, len_seq,bal_seq);

		// complementary node
	bal_word = reverseComplement(word,overlap);
	bal_j = len_seq-0-overlap;  //  0;
	flag = bal_flag = 0;
	if(word<bal_word){
		found = search_kmerset(kset,word,&node);
	}else{
		found = search_kmerset(kset,bal_word,&node);
	}
	if(found&&!node->linear&&!node->checked){
		if(!flag&&node->inEdge==1){
			ctg1start = searchKmerOnCtg(word,KmerCtg1,len1);
			if(ctg1start>0){
				flag = 1;
				seqStart = j + overlap-1;
			}
		}
		if(!bal_flag&&node->inEdge==2){
			bal_ctg2end = searchKmerOnCtg(bal_word,KmerCtg2,len2);
			if(bal_ctg2end>0){
				bal_flag = 2;
				bal_end = bal_j+overlap-1;
			}
		}
	}

	for(j = 1; j <= len_seq - overlap; j ++)	{
		word = nextKmerLocal(word,src_seq[j-1+overlap],WORDF);
		bal_j = len_seq-j-overlap; //  j;
		bal_word = prevKmerLocal(bal_word,bal_seq[bal_j],overlap);

		if(word<bal_word){
			found = search_kmerset(kset,word,&node);
		}else{
			found = search_kmerset(kset,bal_word,&node);
		}
		if(found&&!node->linear&&!node->checked){
			if(!flag&&node->inEdge==1){
				ctg1start = searchKmerOnCtg(word,KmerCtg1,len1);
				if(ctg1start>0){
					flag = 1;
					seqStart = j + overlap-1;
				}
			}else if(flag==1&&node->inEdge==1){
				index = searchKmerOnCtg(word,KmerCtg1,len1);
				if(index>ctg1start){    // choose hit closer to gap
					ctg1start = index;
					seqStart = j + overlap-1;
				}
			}else if(flag==1&&node->inEdge==2){
				ctg2end = searchKmerOnCtg(word,KmerCtg2,len2);
				if(ctg2end>0){
					flag = 3;
					seqEnd = j+overlap-1;
					break;
				}
			}

			if(!bal_flag&&node->inEdge==2){
				bal_ctg2end = searchKmerOnCtg(bal_word,KmerCtg2,len2);
				if(bal_ctg2end>0){
					bal_flag = 2;
					bal_end = bal_j+overlap-1;
				}
			}else if(bal_flag==2&&node->inEdge==2){
				index = searchKmerOnCtg(bal_word,KmerCtg2,len2);
				if(index<bal_ctg2end){   // choose hit closer to gap
					index = bal_ctg2end;
					bal_end = bal_j+overlap-1;
				}
			}else if(bal_flag==2&&node->inEdge==1){
				bal_ctg1start = searchKmerOnCtg(bal_word,KmerCtg1,len1);
				if(bal_ctg1start>0){
					bal_flag = 3;
					bal_start = bal_j+overlap-1;
					break;
				}
			}
		}
	}
	if(flag==3){
		*start = seqStart;
		*end = seqEnd;
		*bal = 0;
		*index1 = ctg1start;
		*index2 = ctg2end;
		return 1;
	}else if(bal_flag==3){
		*start = bal_start;
		*end = bal_end;
		*bal = 1;
		*index1 = bal_ctg1start;
		*index2 = bal_ctg2end;
		return 1;
	}
	return 0;
}

static boolean readsCrossGap(READNEARBY *rdArray, int num, int originOverlap,DARRAY *gapSeqArray,
				Kmer *kmerCtg1,Kmer *kmerCtg2,int overlap,int len1,int len2,
					CTGinSCAF *ctg1,CTGinSCAF *ctg2,KmerSet *kmerS,Kmer WordFilter,int min,int max)
{
	int i,j,start,end,startOnCtg1,endOnCtg2;
	char *bal_seq;
	char *src_seq;
	char *pt;
	boolean bal,ret=0,FILL;

	src_seq = (char *)ckalloc(maxReadLen*sizeof(char));
	bal_seq = (char *)ckalloc(maxReadLen*sizeof(char));

	for(i=0;i<num;i++){
		getSeqFromRead(rdArray[i],src_seq);
		FILL = chopReadFillGap(rdArray[i].len,overlap,src_seq,bal_seq,
					kmerS,WordFilter,&start,&end,&bal,
					kmerCtg1,len1,kmerCtg2,len2,&startOnCtg1,&endOnCtg2);

		if(!FILL||(end-start)<min||(end-start)>max)
			continue;
		fprintf(stderr,"Read across\n");
		//printf("Filled: K %d, ctg1 %d ctg2 %d,start %d end %d\n",overlap,startOnCtg1,endOnCtg2,start,end);
		if(overlap+(len1-startOnCtg1-1)+endOnCtg2-(end-start)>(int)originOverlap)
			continue;  // contig1 and contig2 could not overlap more than origOverlap bases

		ctg1->gapSeqOffset = gapSeqArray->item_c;
		ctg1->gapSeqLen = end-start;
		if(!darrayPut(gapSeqArray,ctg1->gapSeqOffset+(end-start)/4))
			continue;
		pt = (char *)darrayPut(gapSeqArray,ctg1->gapSeqOffset);
		for(j=start+1;j<=end;j++){
			if(bal)
				writeChar2tightString(bal_seq[j],pt,j-start-1);
			else
				writeChar2tightString(src_seq[j],pt,j-start-1);

		}
		ctg1->cutTail = len1-startOnCtg1-1;
		ctg2->cutHead = overlap + endOnCtg2;
		ctg2->scaftig_start = 0;

		ret = 1;
		break;
	}

	free((void*)src_seq);
	free((void*)bal_seq);
	return ret;
}
*/
static void kmerSet_markTandem(KmerSet *set, Kmer WordFilter, int overlap);
static boolean readsCrossGap(READNEARBY *rdArray, int num, int originOverlap, DARRAY *gapSeqArray,
                             Kmer *kmerCtg1, Kmer *kmerCtg2, int overlap,
                             CTGinSCAF *ctg1, CTGinSCAF *ctg2, KmerSet *kmerS, Kmer WordFilter, int min, int max,
                             int offset1, int offset2, char *seqGap, char *seqCtg1, char *seqCtg2, int cut1, int cut2);

int localGraph(READNEARBY *rdArray, int num, CTGinSCAF *ctg1, CTGinSCAF *ctg2,
               int origOverlap, Kmer *kmerCtg1, Kmer *kmerCtg2,
               int overlap, DARRAY *gapSeqArray, char *seqCtg1, char *seqCtg2, char *seqGap)
{
  /**************** put kmer in DBgraph ****************/
  KmerSet *kmerSet;
  Kmer WordFilter = (((Kmer) 1) << (2 * overlap)) - 1;
  /*
  	if(ctg1->ctgID==56410&&ctg2->ctgID==61741)
  		printf("Extract %d reads for gap [%d %d]\n",num,ctg1->ctgID,ctg2->ctgID);
  */
  kmerSet = readsInGap2DBgraph(rdArray, num, ctg1, ctg2, origOverlap,
                               kmerCtg1, kmerCtg2, overlap, WordFilter);
  time_t tt;
  time(&tt);
  //	srand48((int)tt);
  /*
  	int i,j;
  	for(i=0;i<2;i++){
  		int array[] = {0,1,2,3};
  		for(j=4;j>0;j--)
  			fprintf(stderr,"%d ", nPick1(array,j));
  	}
  	fprintf(stderr,"\n");
  */
  /***************** search path to connect contig ends ********/
  int gapLen = ctg2->start - ctg1->end - origOverlap + overlap;
  int min = gapLen - GLDiff > 0 ? gapLen - GLDiff : 0;
  int max = gapLen + GLDiff < 10 ? 10 : gapLen + GLDiff;
  //count kmer number for contig1 and contig2 ends
  int len1, len2;
  len1 = CTGendLen < contig_array[ctg1->ctgID].length + origOverlap ?
         CTGendLen : contig_array[ctg1->ctgID].length + origOverlap;
  len2 = CTGendLen < contig_array[ctg2->ctgID].length + origOverlap ?
         CTGendLen : contig_array[ctg2->ctgID].length + origOverlap;
  len1 -= overlap - 1;
  len2 -= overlap - 1;

  //int pathNum = 2;
  int offset1 = 0, offset2 = 0, cut1 = 0, cut2 = 0;
  int pathNum = searchFgap(kmerSet, ctg1, ctg2, kmerCtg1, kmerCtg2,
                           origOverlap, overlap, gapSeqArray,
                           len1, len2, WordFilter, &offset1, &offset2, seqGap, &cut1, &cut2);

  //printf("SF: %d K %d\n",pathNum,overlap);
  if(pathNum == 0)
    {
      free_kmerset(kmerSet);
      return 0;
    }
  else if(pathNum == 1)
    {
      free_kmerset(kmerSet);
      return 1;
    }/*

	else{
		printf("ret %d\n",pathNum);
		free_kmerset(kmerSet);
		return 0;
	}  */

  /******************* cross the gap by single reads *********/
  //kmerSet_markTandem(kmerSet,WordFilter,overlap);
  maskRepeatNode(kmerSet, kmerCtg1, kmerCtg2, overlap,
                 len1, len2, max, WordFilter);
  boolean found = readsCrossGap(rdArray, num, origOverlap, gapSeqArray,
                                kmerCtg1, kmerCtg2, overlap, ctg1, ctg2, kmerSet, WordFilter, min, max,
                                offset1, offset2, seqGap, seqCtg1, seqCtg2, cut1, cut2);

  if(found)
    {
      //fprintf(stderr,"read across\n");
      free_kmerset(kmerSet);
      return found;
    }
  else
    {
      free_kmerset(kmerSet);
      return 0;
    }

}

static void kmerSet_markTandem(KmerSet *set, Kmer WordFilter, int overlap)
{
  kmer_t *rs;
  long long counter = 0;
  int num_route, steps;
  int min = 1, max = overlap, maxRoute = 1;
  int traceCounter;

  set->iter_ptr = 0;

  while(set->iter_ptr < set->size)
    {
      if(!is_kmer_entity_null(set->flags, set->iter_ptr))
        {
          rs = set->array + set->iter_ptr;

          if(rs->inEdge > 0)
            {
              set->iter_ptr ++;
              continue;
            }

          num_route = traceCounter = 0;
          steps = 0;
          trace4Repeat(rs->seq, steps, min, max, &num_route,
                       set, rs->seq, overlap, WordFilter,
                       &traceCounter, maxRoute, NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL);

          if(num_route < 1)
            {
              set->iter_ptr ++;
              continue;
            }

          /*
          printKmer(stderr,rs->seq,overlap);
          fprintf(stderr, "\n");
          */
          rs->checked = 1;
          counter++;
        }

      set->iter_ptr ++;
    }
}
/******************* the following is for read-crossing gaps *************************/

#define MAXREADLENGTH 100

static const int INDEL = 0;
static const int SIM[4][4] =
{
  {1, 0, 0, 0},
  {0, 1, 0, 0},
  {0, 0, 1, 0},
  {0, 0, 0, 1}
};
static char fastSequence[MAXREADLENGTH];
static char slowSequence[MAXREADLENGTH];

static int Fmatrix[MAXREADLENGTH + 1][MAXREADLENGTH + 1];
static int slowToFastMapping[MAXREADLENGTH + 1];
static int fastToSlowMapping[MAXREADLENGTH + 1];

static int max(int A, int B, int C)
{
  A = A >= B ? A : B;
  return (A >= C ? A : C);

}

static int compareSequences(char *sequence1, char *sequence2, int length1, int length2)
{
  if(length1 < 1 || length2 < 1 || length1 > MAXREADLENGTH || length2 > MAXREADLENGTH)
    return 0;

  int i, j;
  int Choice1, Choice2, Choice3;
  int maxScore;

  for (i = 0; i <= length1; i++)
    Fmatrix[i][0] = 0;

  for (j = 0; j <= length2; j++)
    Fmatrix[0][j] = 0;

  for (i = 1; i <= length1; i++)
    {
      for (j = 1; j <= length2; j++)
        {
          Choice1 =
            Fmatrix[i - 1][j - 1] +
            SIM[(int) sequence1[i - 1]]
            [(int) sequence2[j - 1]];
          Choice2 = Fmatrix[i - 1][j] + INDEL;
          Choice3 = Fmatrix[i][j - 1] + INDEL;
          Fmatrix[i][j] = max(Choice1, Choice2, Choice3);
        }
    }

  maxScore = Fmatrix[length1][length2];
  return maxScore;
}

static void mapSlowOntoFast(int slowSeqLength, int fastSeqLength)
{
  int slowIndex = slowSeqLength;
  int fastIndex = fastSeqLength;
  int fastn, slown;

  if (slowIndex == 0)
    {
      slowToFastMapping[0] = fastIndex;

      while (fastIndex >= 0)
        fastToSlowMapping[fastIndex--] = 0;

      return;
    }

  if (fastIndex == 0)
    {
      while (slowIndex >= 0)
        slowToFastMapping[slowIndex--] = 0;

      fastToSlowMapping[0] = slowIndex;

      return;
    }

  while (slowIndex > 0 && fastIndex > 0)
    {
      fastn = (int) fastSequence[fastIndex - 1]; //getCharInTightString(fastSequence,fastIndex-1);
      slown = (int) slowSequence[slowIndex - 1]; //getCharInTightString(slowSequence,slowIndex-1);

      if (Fmatrix[fastIndex][slowIndex] ==
          Fmatrix[fastIndex - 1][slowIndex - 1] +
          SIM[fastn][slown])
        {
          fastToSlowMapping[--fastIndex] = --slowIndex;
          slowToFastMapping[slowIndex] = fastIndex;
        }
      else if (Fmatrix[fastIndex][slowIndex] ==
               Fmatrix[fastIndex - 1][slowIndex] + INDEL)
        fastToSlowMapping[--fastIndex] = slowIndex - 1;

      else if (Fmatrix[fastIndex][slowIndex] ==
               Fmatrix[fastIndex][slowIndex - 1] + INDEL)
        slowToFastMapping[--slowIndex] = fastIndex - 1;

      else
        {
          printf("compareSequence: Error trace\n");
          fflush(stdout);
          abort();
        }
    }

  while (slowIndex > 0)
    slowToFastMapping[--slowIndex] = -1;

  while (fastIndex > 0)
    fastToSlowMapping[--fastIndex] = -1;

  slowToFastMapping[slowSeqLength] =
    fastSeqLength;
  fastToSlowMapping[fastSeqLength] =
    slowSeqLength;
}

static boolean chopReadFillGap(int len_seq, int overlap, char *src_seq, char *bal_seq,
                               KmerSet *kset, Kmer WORDF, int *start, int *end, boolean *bal,
                               Kmer *KmerCtg1, int len1, Kmer *KmerCtg2, int len2, int *index1, int *index2)
{
  int index, j = 0, bal_j;
  Kmer word, bal_word;
  int flag = 0, bal_flag = 0;
  int ctg1start, bal_ctg1start, ctg2end, bal_ctg2end;
  int seqStart, bal_start, seqEnd, bal_end;
  kmer_t *node;
  boolean found;

  if(len_seq < overlap + 1)
    {
      return 0;
    }

  word = 0;

  for (index = 0; index < overlap; index++)
    {
      word <<= 2;
      word += src_seq[index];
    }

  reverseComplementSeq(src_seq, len_seq, bal_seq);

  // complementary node
  bal_word = reverseComplement(word, overlap);
  bal_j = len_seq - 0 - overlap; //  0;
  flag = bal_flag = 0;

  if(word < bal_word)
    {
      found = search_kmerset(kset, word, &node);
    }
  else
    {
      found = search_kmerset(kset, bal_word, &node);
    }

  if(found && !node->linear && !node->checked)
    {
      if(!flag && node->inEdge == 1)
        {
          ctg1start = searchKmerOnCtg(word, KmerCtg1, len1);

          if(ctg1start >= 0)
            {
              flag = 1;
              seqStart = j + overlap - 1;
            }
        }

      if(!bal_flag && node->inEdge == 2)
        {
          bal_ctg2end = searchKmerOnCtg(bal_word, KmerCtg2, len2);

          if(bal_ctg2end >= 0)
            {
              bal_flag = 2;
              bal_end = bal_j + overlap - 1;
            }
        }
    }

  for(j = 1; j <= len_seq - overlap; j ++)
    {
      word = nextKmerLocal(word, src_seq[j - 1 + overlap], WORDF);
      bal_j = len_seq - j - overlap; //  j;
      bal_word = prevKmerLocal(bal_word, bal_seq[bal_j], overlap);

      if(word < bal_word)
        {
          found = search_kmerset(kset, word, &node);
        }
      else
        {
          found = search_kmerset(kset, bal_word, &node);
        }

      if(found && !node->linear && !node->checked)
        {
          if(!flag && node->inEdge == 1)
            {
              ctg1start = searchKmerOnCtg(word, KmerCtg1, len1);

              if(ctg1start >= 0)
                {
                  flag = 1;
                  seqStart = j + overlap - 1;
                }
            }
          else if(flag == 1 && node->inEdge == 1)
            {
              index = searchKmerOnCtg(word, KmerCtg1, len1);

              if(index >= 0 && index > ctg1start) // choose hit closer to gap
                {
                  ctg1start = index;
                  seqStart = j + overlap - 1;
                }
            }
          else if(flag == 1 && node->inEdge == 2)
            {
              ctg2end = searchKmerOnCtg(word, KmerCtg2, len2);

              if(ctg2end >= 0)
                {
                  flag = 3;
                  seqEnd = j + overlap - 1;
                  break;
                }
            }

          if(!bal_flag && node->inEdge == 2)
            {
              bal_ctg2end = searchKmerOnCtg(bal_word, KmerCtg2, len2);

              if(bal_ctg2end >= 0)
                {
                  bal_flag = 2;
                  bal_end = bal_j + overlap - 1;
                }
            }
          else if(bal_flag == 2 && node->inEdge == 2)
            {
              index = searchKmerOnCtg(bal_word, KmerCtg2, len2);

              if(index >= 0 && index < bal_ctg2end) // choose hit closer to gap
                {
                  bal_ctg2end = index;
                  bal_end = bal_j + overlap - 1;
                }
            }
          else if(bal_flag == 2 && node->inEdge == 1)
            {
              bal_ctg1start = searchKmerOnCtg(bal_word, KmerCtg1, len1);

              if(bal_ctg1start >= 0)
                {
                  bal_flag = 3;
                  bal_start = bal_j + overlap - 1;
                  break;
                }
            }
        }
    }

  if(flag == 3)
    {
      *start = seqStart;
      *end = seqEnd;
      *bal = 0;
      *index1 = ctg1start;
      *index2 = ctg2end;
      return 1;
    }
  else if(bal_flag == 3)
    {
      *start = bal_start;
      *end = bal_end;
      *bal = 1;
      *index1 = bal_ctg1start;
      *index2 = bal_ctg2end;
      return 1;
    }

  return 0;
}


static int cutSeqFromTightStr(char *tightStr, int length, int start, int end, int revS, char *src_seq)
{
  int i, index = 0;
  end = end < length ? end : length - 1;
  start = start >= 0 ? start : 0;

  if(!revS)
    {
      for(i = start; i <= end; i++)
        src_seq[index++] = getCharInTightString(tightStr, i);
    }
  else
    {
      for(i = length - 1 - start; i >= length - end - 1; i--)
        src_seq[index++] = int_comp(getCharInTightString(tightStr, i));
    }

  return end - start + 1;
}

static int cutSeqFromCtg(unsigned int ctgID, int start, int end, char *sequence, int originOverlap)
{

  unsigned int bal_ctg = getTwinCtg(ctgID);

  if(contig_array[ctgID].length < 1)
    return 0;

  int length = contig_array[ctgID].length + originOverlap;

  if(contig_array[ctgID].seq)
    return cutSeqFromTightStr(contig_array[ctgID].seq, length, start, end, 0, sequence);
  else
    return cutSeqFromTightStr(contig_array[bal_ctg].seq, length, start, end, 1, sequence);

}

static int cutSeqFromRead(char *src_seq, int length, int start, int end, char *sequence)
{
  if(end >= length)
    printf("******: end %d length %d\n", end, length);

  end = end < length ? end : length - 1;
  start = start >= 0 ? start : 0;
  int i;

  for(i = start; i <= end; i++)
    sequence[i - start] = src_seq[i];

  return end - start + 1;
}

static void printSeq(FILE *fo, char *seq, int len)
{
  int i;

  for(i = 0; i < len; i++)
    fprintf(fo, "%c", int2base((int)seq[i]));

  fprintf(fo, "\n");
}

static boolean readsCrossGap(READNEARBY *rdArray, int num, int originOverlap, DARRAY *gapSeqArray,
                             Kmer *kmerCtg1, Kmer *kmerCtg2, int overlap,
                             CTGinSCAF *ctg1, CTGinSCAF *ctg2, KmerSet *kmerS, Kmer WordFilter, int min, int max,
                             int offset1, int offset2, char *seqGap, char *seqCtg1, char *seqCtg2, int cut1, int cut2)
{
  int i, j, start, end, startOnCtg1, endOnCtg2;
  char *bal_seq;
  char *src_seq;
  char *pt;
  boolean bal, ret = 0, FILL;
  double maxScore = 0.0;
  int maxIndex;
  int lenCtg1, lenCtg2;
  //build sequences on left and right of the uncertain region
  int buffer_size = maxReadLen > 100 ? maxReadLen : 100;
  int length = contig_array[ctg1->ctgID].length + originOverlap;

  if(buffer_size > offset1)
    {
      lenCtg1 = cutSeqFromCtg(ctg1->ctgID, length - cut1 - (buffer_size - offset1), length - 1 - cut1, seqCtg1, originOverlap);

      for(i = 0; i < offset1; i++)
        seqCtg1[lenCtg1 + i] = seqGap[i];

      lenCtg1 += offset1;
    }
  else
    {
      for(i = offset1 - buffer_size; i < offset1; i++)
        seqCtg1[i + buffer_size - offset1] = seqGap[i];

      lenCtg1 = buffer_size;
    }

  length = contig_array[ctg2->ctgID].length + originOverlap;

  if(buffer_size > offset2)
    {
      lenCtg2 = cutSeqFromCtg(ctg2->ctgID, cut2, buffer_size - offset2 - 1 + cut2, &(seqCtg2[offset2]), originOverlap);

      for(i = 0; i < offset2; i++)
        seqCtg2[i] = seqGap[i + offset1];

      lenCtg2 += offset2;
    }
  else
    {
      for(i = 0; i < buffer_size; i++)
        seqCtg2[i] = seqGap[i + offset1];

      lenCtg2 = buffer_size;
    }

  /*
  if(offset1>0||offset2>0){
  for(i=0;i<lenCtg1;i++)
  	fprintf(stderr,"%c",int2base(seqCtg1[i]));
  fprintf(stderr,": CTG1\n");
  for(i=0;i<lenCtg2;i++)
  	fprintf(stderr,"%c",int2base(seqCtg2[i]));
  fprintf(stderr,": CTG2\n");
  }
  */
  //chop kmer from both ends of the uncertain region
  int len1, len2;
  len1 = CTGendLen < lenCtg1 ? CTGendLen : lenCtg1;
  len2 = CTGendLen < lenCtg2 ? CTGendLen : lenCtg2;
  chopKmer4Ctg(kmerCtg1, len1, overlap, &(seqCtg1[lenCtg1 - len1]), WordFilter);
  chopKmer4Ctg(kmerCtg2, len2, overlap, seqCtg2, WordFilter);
  len1 -= overlap - 1;
  len2 -= overlap - 1;

  src_seq = (char *)ckalloc(maxReadLen * sizeof(char));
  bal_seq = (char *)ckalloc(maxReadLen * sizeof(char));

  int *START = (int *)ckalloc(num * sizeof(int));
  int *END = (int *)ckalloc(num * sizeof(int));
  int *INDEX1 = (int *)ckalloc(num * sizeof(int));
  int *INDEX2 = (int *)ckalloc(num * sizeof(int));
  double *SCORE = (double *)ckalloc(num * sizeof(double));
  boolean *BAL = (boolean *)ckalloc(num * sizeof(boolean));
  memset(SCORE, 0, num * sizeof(double));

  for(i = 0; i < num; i++)
    {
      getSeqFromRead(rdArray[i], src_seq);
      FILL = chopReadFillGap(rdArray[i].len, overlap, src_seq, bal_seq,
                             kmerS, WordFilter, &start, &end, &bal,
                             kmerCtg1, len1, kmerCtg2, len2, &startOnCtg1, &endOnCtg2);

      if(!FILL || (end - start) < min || (end - start) > max)
        continue;

      if(overlap + (len1 - startOnCtg1 - 1) + endOnCtg2 - (end - start) > (int)originOverlap)
        continue;  // contig1 and contig2 could not overlap more than origOverlap bases

      START[i] = start;
      END[i] = end;
      INDEX1[i] = startOnCtg1;
      INDEX2[i] = endOnCtg2;
      BAL[i] = bal;

      int matchLen = 2 * overlap < (end - start + overlap) ? 2 * overlap : (end - start + overlap);
      int match;
      int alignLen = matchLen;
      //compare the left of hit kmer on ctg1
      //int ctgLeft = (contig_array[ctg1->ctgID].length+originOverlap)-(len1+overlap-1)+startOnCtg1;
      int ctgLeft = (lenCtg1) - (len1 + overlap - 1) + startOnCtg1;
      int readLeft = start - overlap + 1;
      int cmpLen = ctgLeft < readLeft ? ctgLeft : readLeft;

      cmpLen = cmpLen <= MAXREADLENGTH ? cmpLen : MAXREADLENGTH;
      //cutSeqFromCtg(ctg1->ctgID,ctgLeft-cmpLen,ctgLeft-1,fastSequence,originOverlap);
      cutSeqFromRead(seqCtg1, lenCtg1, ctgLeft - cmpLen, ctgLeft - 1, fastSequence);

      if(!bal)
        cutSeqFromRead(src_seq, rdArray[i].len, readLeft - cmpLen, readLeft - 1, slowSequence);
      else
        cutSeqFromRead(bal_seq, rdArray[i].len, readLeft - cmpLen, readLeft - 1, slowSequence);

      match = compareSequences(fastSequence, slowSequence, cmpLen, cmpLen);

      alignLen += cmpLen;
      matchLen += match;

      //compare the right of hit kmer on ctg1
      int ctgRight = len1 - startOnCtg1 - 1;

      cmpLen = ctgRight < (rdArray[i].len - start - 1) ? ctgRight : (rdArray[i].len - start - 1);
      cmpLen = cmpLen <= MAXREADLENGTH ? cmpLen : MAXREADLENGTH;
      //cutSeqFromCtg(ctg1->ctgID,ctgLeft+overlap,ctgLeft+overlap+cmpLen-1,fastSequence,originOverlap);
      cutSeqFromRead(seqCtg1, lenCtg1, ctgLeft + overlap, ctgLeft + overlap + cmpLen - 1, fastSequence);

      if(!bal)
        cutSeqFromRead(src_seq, rdArray[i].len, start + 1, start + cmpLen, slowSequence);
      else
        cutSeqFromRead(bal_seq, rdArray[i].len, start + 1, start + cmpLen, slowSequence);

      match = compareSequences(fastSequence, slowSequence, cmpLen, cmpLen);
      //fprintf(stderr,"%d -- %d\n",match,cmpLen);

      alignLen += cmpLen;
      matchLen += match;

      //compare the left of hit kmer on ctg2
      ctgLeft = endOnCtg2;
      readLeft = end - overlap + 1;
      cmpLen = ctgLeft < readLeft ? ctgLeft : readLeft;
      cmpLen = ctgLeft <= MAXREADLENGTH ? ctgLeft : MAXREADLENGTH;
      //cutSeqFromCtg(ctg2->ctgID,endOnCtg2-cmpLen,endOnCtg2-1,fastSequence,originOverlap);
      cutSeqFromRead(seqCtg2, lenCtg2, endOnCtg2 - cmpLen, endOnCtg2 - 1, fastSequence);

      if(!bal)
        cutSeqFromRead(src_seq, rdArray[i].len, readLeft - cmpLen, readLeft - 1, slowSequence);
      else
        cutSeqFromRead(bal_seq, rdArray[i].len, readLeft - cmpLen, readLeft - 1, slowSequence);

      match = compareSequences(fastSequence, slowSequence, cmpLen, cmpLen);
      alignLen += cmpLen;
      matchLen += match;

      //compare the right of hit kmer on ctg2
      //ctgRight = contig_array[ctg2->ctgID].length+originOverlap-endOnCtg2-overlap;
      ctgRight = lenCtg2 - endOnCtg2 - overlap;
      cmpLen = ctgRight < (rdArray[i].len - end - 1) ? ctgRight : (rdArray[i].len - end - 1);
      cmpLen = cmpLen <= MAXREADLENGTH ? cmpLen : MAXREADLENGTH;
      //cutSeqFromCtg(ctg2->ctgID,endOnCtg2+overlap,endOnCtg2+overlap+cmpLen-1,fastSequence,originOverlap);
      cutSeqFromRead(seqCtg2, lenCtg2, endOnCtg2 + overlap, endOnCtg2 + overlap + cmpLen - 1, fastSequence);

      if(!bal)
        cutSeqFromRead(src_seq, rdArray[i].len, end + 1, end + cmpLen, slowSequence);
      else
        cutSeqFromRead(bal_seq, rdArray[i].len, end + 1, end + cmpLen, slowSequence);

      match = compareSequences(fastSequence, slowSequence, cmpLen, cmpLen);
      alignLen += cmpLen;
      matchLen += match;
      /*
      if(cmpLen>0&&match!=cmpLen+overlap){
      	printSeq(stderr,fastSequence,cmpLen+overlap);
      	printSeq(stderr,slowSequence,cmpLen+overlap);
      	printKmer(stderr,kmerCtg2[endOnCtg2],overlap);
      	fprintf(stderr,": %d(%d)\n",bal,endOnCtg2);
      }else if(cmpLen>0&&match==cmpLen+overlap)
      	fprintf(stderr,"Perfect\n");
      */
      double score = (double)matchLen / alignLen;

      if(maxScore < score)
        {
          maxScore = score;
          //fprintf(stderr,"%4.2f (%d/%d)\n",maxScore,matchLen,alignLen);
          maxIndex = i;
        }

      SCORE[i] = score;
    }

  /*
  if(maxScore>0.0)
  	fprintf(stderr,"SCORE: %4.2f\n",maxScore);
  */
  if(maxScore > 0.9)
    {
      /*
      for(i=0;i<lenCtg1;i++)
      	fprintf(stderr,"%c",int2base(seqCtg1[i]));
      fprintf(stderr,": CTG1\n");
      for(i=0;i<lenCtg2;i++)
      	fprintf(stderr,"%c",int2base(seqCtg2[i]));
      fprintf(stderr,": CTG2\n");
      fprintf(stderr,"%d+%d -- %d+%d, SCORE: %4.2f\n ",offset1,offset2,cut1,cut2,maxScore);
      */
      getSeqFromRead(rdArray[maxIndex], src_seq);
      reverseComplementSeq(src_seq, rdArray[maxIndex].len, bal_seq);

      int leftRemain = offset1 - (len1 - INDEX1[maxIndex] - 1) > 0 ? offset1 - (len1 - INDEX1[maxIndex] - 1) : 0;
      int rightRemain = offset2 - (overlap + INDEX2[maxIndex]) > 0 ? offset2 - (overlap + INDEX2[maxIndex]) : 0;

      ctg1->gapSeqOffset = gapSeqArray->item_c;
      ctg1->gapSeqLen = END[maxIndex] - START[maxIndex] + leftRemain + rightRemain;

      if(darrayPut(gapSeqArray, ctg1->gapSeqOffset + (END[maxIndex] - START[maxIndex] + leftRemain + rightRemain) / 4))
        {
          pt = (char *)darrayPut(gapSeqArray, ctg1->gapSeqOffset);

          for(j = 0; j < leftRemain; j++) //get the left side of the gap region from search
            {
              writeChar2tightString(seqGap[j], pt, j);
              fprintf(stderr, "%c", int2base(seqGap[j]));
            }

          for(j = START[maxIndex] + 1; j <= END[maxIndex]; j++)
            {
              if(BAL[maxIndex])
                {
                  writeChar2tightString(bal_seq[j], pt, j - START[maxIndex] - 1 + leftRemain);
                  fprintf(stderr, "%c", int2base(bal_seq[j]));
                }
              else
                {
                  writeChar2tightString(src_seq[j], pt, j - START[maxIndex] - 1 + leftRemain);
                  fprintf(stderr, "%c", int2base(src_seq[j]));
                }
            }

          for(j = offset2 - rightRemain; j < offset2; j++) //get the right side of the gap region from search
            {
              writeChar2tightString(seqGap[j + leftRemain], pt, j + END[maxIndex] - START[maxIndex] + leftRemain);
              fprintf(stderr, "%c", int2base(seqGap[j + leftRemain]));
            }

          fprintf(stderr, ": GAPSEQ (%d+%d)(%d+%d)(%d+%d)(%d+%d) B %d\n", offset1, offset2, cut1, cut2,
                  len1 - INDEX1[maxIndex] - 1, INDEX2[maxIndex], START[maxIndex], END[maxIndex], BAL[maxIndex]);

          ctg1->cutTail = len1 - INDEX1[maxIndex] - 1 - offset1 + cut1 > cut1 ? len1 - INDEX1[maxIndex] - 1 - offset1 + cut1 : cut1;
          ctg2->cutHead = overlap + INDEX2[maxIndex] - offset2 + cut2 > cut2 ? overlap + INDEX2[maxIndex] - offset2 + cut2 : cut2;
          ctg2->scaftig_start = 0;
          ret = 1;
        }
    }

  free((void *)START);
  free((void *)END);
  free((void *)INDEX1);
  free((void *)INDEX2);
  free((void *)SCORE);
  free((void *)BAL);

  free((void *)src_seq);
  free((void *)bal_seq);
  return ret;
}