static char rcsid[] = "$Id: indel.c 204387 2017-03-18 00:02:54Z twu $";
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "indel.h"

#include "assert.h"
#include "mem.h"
#include "genome128_hr.h"
#include "stage3hr.h"
#include "intron.h"


/* Indels */ 
#ifdef DEBUG2
#define debug2(x) x
#else
#define debug2(x)
#endif


static int min_indel_end_matches;
static int indel_penalty_middle;


void
Indel_setup (int min_indel_end_matches_in, int indel_penalty_middle_in) {
  min_indel_end_matches = min_indel_end_matches_in;
  indel_penalty_middle = indel_penalty_middle_in;
  return;
}


/* Called only by sarray-read.c, where plusp is always true */
/* indels is positive here */
int
Indel_resolve_middle_insertion (int *best_nmismatches_i, int *best_nmismatches_j,
				Univcoord_T left, int indels, Compress_T query_compress,
				int querystart, int queryend, int querylength,
				int max_mismatches_allowed, bool plusp, int genestrand) {
  int best_indel_pos = -1, indel_pos;
#ifdef DEBUG2
  int i;
  char *gbuffer;
#endif
  int nmismatches_left, nmismatches_right;
  int best_sum, sum, nmismatches_lefti, nmismatches_righti, lefti, righti;
  int *mismatch_positions_left, *mismatch_positions_right;

#ifdef HAVE_ALLOCA
  if (querylength <= MAX_STACK_READLENGTH) {
    mismatch_positions_left = (int *) ALLOCA(querylength * sizeof(int));
    mismatch_positions_right = (int *) ALLOCA(querylength * sizeof(int));
  } else {
    mismatch_positions_left = (int *) MALLOC(querylength * sizeof(int));
    mismatch_positions_right = (int *) MALLOC(querylength * sizeof(int));
  }
#else
  mismatch_positions_left = (int *) MALLOC(querylength * sizeof(int));
  mismatch_positions_right = (int *) MALLOC(querylength * sizeof(int));
#endif

  if (max_mismatches_allowed > querylength) {
    max_mismatches_allowed = querylength;
  }

  /* query has insertion.  Get |indels| less from genome; trim from left. */
  /* left = ptr->diagonal - querylength; */

  assert(indels > 0);
  debug2(gbuffer = (char *) CALLOC(querylength-indels+1,sizeof(char)));
  debug2(Genome_fill_buffer_blocks(left+indels,querylength-indels,gbuffer));
  debug2(printf("solve_middle_indel, plus, insertion: Getting genome at diagonal - querylength %d + indels %d = %llu\n",
		querylength,indels,(unsigned long long) left+indels));
  debug2(printf("g1: %s\n",gbuffer));
  debug2(printf("g2: %s\n",&(gbuffer[indels])));

  /* No need to check chromosome bounds */
  debug2(printf("max_mismatches_allowed is %d.  Calling Genome_mismatches_left over %d..%d\n",
		max_mismatches_allowed,querystart,queryend));
  nmismatches_left = Genome_mismatches_left(mismatch_positions_left,max_mismatches_allowed,
					    query_compress,left,/*pos5*/querystart,/*pos3*/queryend,
					    plusp,genestrand);

  debug2(
	 printf("%d mismatches on left at:",nmismatches_left);
	 for (i = 0; i <= nmismatches_left; i++) {
	   printf(" %d",mismatch_positions_left[i]);
	 }
	 printf("\n");
	 );


  /* No need to check chromosome bounds */
  debug2(printf("max_mismatches_allowed is %d.  Calling Genome_mismatches_right over %d..%d\n",
		max_mismatches_allowed,querystart,queryend));
  nmismatches_right = Genome_mismatches_right(mismatch_positions_right,max_mismatches_allowed,
					      query_compress,left-indels,/*pos5*/querystart,/*pos3*/queryend,
					      plusp,genestrand);

  debug2(
	 printf("%d mismatches on right at:",nmismatches_right);
	 for (i = 0; i <= nmismatches_right; i++) {
	   printf(" %d",mismatch_positions_right[i]);
	 }
	 printf("\n");
	 );

  best_sum = querylength + querylength;

  /* Modeled after end D to get lowest possible coordinate */
  righti = 0;
  lefti = nmismatches_left - 1;
  nmismatches_righti = /*righti*/ 0;
  nmismatches_lefti = /*lefti+1*/ nmismatches_left;

  while (righti < nmismatches_right) {
    while (lefti >= 0 && mismatch_positions_left[lefti] > mismatch_positions_right[righti] - indels) {
      lefti--;
    }
    sum = righti + lefti + 1;
    debug2(printf("  (Case D) sum %d=%d+%d at indel_pos %d.",
		  sum,righti,lefti+1,mismatch_positions_right[righti]-indels+1));
    if (sum <= best_sum) {
      indel_pos = mismatch_positions_right[righti] - indels + 1;
      if (indel_pos >= min_indel_end_matches && indel_pos + indels <= querylength - min_indel_end_matches &&
	  /* Require more matches on both ends than indel length */ indel_pos >= indels && indel_pos + indels <= querylength - indels) {
	best_indel_pos = indel_pos;
	nmismatches_righti = righti;
	nmismatches_lefti = lefti + 1;
	debug2(printf("**"));
	best_sum = sum;
      }
    }
    righti++;
  }
  debug2(printf("\n"));


  /* Try from other side to see if we missed anything */
  lefti = 0;
  righti = nmismatches_right - 1;

  while (lefti < nmismatches_left) {
    while (righti >= 0 && mismatch_positions_right[righti] < mismatch_positions_left[lefti] + indels) {
      righti--;
    }
    sum = lefti + righti + 1;
    debug2(printf("  (Case D2) sum %d=%d+%d at indel_pos %d.",
		  sum,lefti,righti+1,mismatch_positions_left[lefti]));
    if (sum < best_sum) {
      indel_pos = mismatch_positions_left[lefti];
      if (indel_pos >= min_indel_end_matches && indel_pos + indels <= querylength - min_indel_end_matches &&
	  /* Require more matches on both ends than insertion length */ indel_pos >= indels && indel_pos + indels <= querylength - indels) {
	best_indel_pos = indel_pos;
	nmismatches_righti = righti + 1;
	nmismatches_lefti = lefti;
	debug2(printf("**"));
	best_sum = sum;
      }
    } else if (sum == best_sum) {
      indel_pos = mismatch_positions_left[lefti];
      if (indel_pos < best_indel_pos) {
	if (indel_pos >= min_indel_end_matches && indel_pos + indels <= querylength - min_indel_end_matches &&
	  /* Require more matches on both ends than insertion length */ indel_pos >= indels && indel_pos + indels <= querylength - indels) {
	  best_indel_pos = indel_pos;
	  nmismatches_righti = righti + 1;
	  nmismatches_lefti = lefti;
	  debug2(printf("**"));
	  /* best_sum = sum; */
	}
      }
    }
    lefti++;
  }
  debug2(printf("\n"));

#ifdef HAVE_ALLOCA
  if (querylength <= MAX_STACK_READLENGTH) {
    FREEA(mismatch_positions_left);
    FREEA(mismatch_positions_right);
  } else {
    FREE(mismatch_positions_left);
    FREE(mismatch_positions_right);
  }
#else
  FREE(mismatch_positions_left);
  FREE(mismatch_positions_right);
#endif

  *best_nmismatches_i = nmismatches_lefti;
  *best_nmismatches_j = nmismatches_righti;

  if (best_sum > max_mismatches_allowed) {
    debug2(printf("Returning -1\n"));
    return -1;
#if 0
  } else if (plusp == true) {
    return best_indel_pos;
  } else {
    return querylength - best_indel_pos - indels;
#else
  } else {
    debug2(printf("Returning %d with mismatches %d+%d\n",
		  best_indel_pos,*best_nmismatches_i,*best_nmismatches_j));
    return best_indel_pos;
#endif
  }
}


/* Called only by sarray-read.c, where plusp is always true */
/* indels is negative here */
int
Indel_resolve_middle_deletion (int *best_introntype,
			       int *best_nmismatches_i, int *best_nmismatches_j,
			       Univcoord_T left, int indels, Compress_T query_compress,
			       int querystart, int queryend, int querylength,
			       int max_mismatches_allowed, bool plusp, int genestrand,
			       int min_intronlength) {
  int best_indel_pos = -1, indel_pos;
  char *gbuffer;
#ifdef DEBUG2
  int i;
#endif
  int nmismatches_left, nmismatches_right, nmismatches_lefti, nmismatches_righti;
  int best_sum, sum, lefti, righti;
  int *mismatch_positions_left, *mismatch_positions_right;
  char left1, left2, right2, right1;
  int introntype, intron_level, best_intron_level;

  *best_introntype = NONINTRON;
  best_intron_level = Intron_level(NONINTRON);

#ifdef HAVE_ALLOCA
  if (querylength <= MAX_STACK_READLENGTH) {
    mismatch_positions_left = (int *) ALLOCA(querylength * sizeof(int));
    mismatch_positions_right = (int *) ALLOCA(querylength * sizeof(int));
  } else {
    mismatch_positions_left = (int *) MALLOC(querylength * sizeof(int));
    mismatch_positions_right = (int *) MALLOC(querylength * sizeof(int));
  }
#else
  mismatch_positions_left = (int *) MALLOC(querylength * sizeof(int));
  mismatch_positions_right = (int *) MALLOC(querylength * sizeof(int));
#endif

  if (max_mismatches_allowed > querylength) {
    max_mismatches_allowed = querylength;
  }

  /* query has deletion.  Get |indels| more from genome; add to right. */
  /* left = ptr->diagonal - querylength; */

  assert(indels < 0);
#ifdef DEBUG2
  gbuffer = (char *) CALLOC(querylength-indels+1,sizeof(char));
  Genome_fill_buffer_blocks(left,querylength-indels,gbuffer);
#else
  if (-indels >= min_intronlength) {
    gbuffer = (char *) CALLOC(querylength-indels+1,sizeof(char));
    Genome_fill_buffer_blocks(left,querylength-indels,gbuffer);
  }  
#endif
  debug2(printf("solve_middle_indel, plus, deletion (indels %d), max_mismatches_allowed %d: Getting genome at diagonal - querylength %d = %llu\n",
		indels,max_mismatches_allowed,querylength,(unsigned long long) left));
  debug2(printf("g1: %s\n",gbuffer));
  debug2(printf("g2: %s\n",&(gbuffer[-indels])));

  /* No need to check chromosome bounds */
  nmismatches_left = Genome_mismatches_left(mismatch_positions_left,max_mismatches_allowed,
					    query_compress,left,/*pos5*/querystart,/*pos3*/queryend,
					    plusp,genestrand);

  debug2(
	 printf("%d mismatches on left at:",nmismatches_left);
	 for (i = 0; i <= nmismatches_left; i++) {
	   printf(" %d",mismatch_positions_left[i]);
	 }
	 printf("\n");
	 );

  /* No need to check chromosome bounds */
  nmismatches_right = Genome_mismatches_right(mismatch_positions_right,max_mismatches_allowed,
					      query_compress,left-indels,/*pos5*/querystart,/*pos3*/queryend,
					      plusp,genestrand);

  debug2(
	 printf("%d mismatches on right at:",nmismatches_right);
	 for (i = 0; i <= nmismatches_right; i++) {
	   printf(" %d",mismatch_positions_right[i]);
	 }
	 printf("\n");
	 );

  best_sum = querylength + querylength;

  /* Modeled after end C to get lowest possible coordinate */
  righti = 0;
  lefti = nmismatches_left - 1;
  nmismatches_righti = /*righti*/ 0;
  nmismatches_lefti = /*lefti+1*/ nmismatches_left;

  while (righti < nmismatches_right) {
    while (lefti >= 0 && mismatch_positions_left[lefti] > mismatch_positions_right[righti]) {
      lefti--;
    }
    sum = righti + lefti + 1;

    if (-indels >= min_intronlength) {
      /* Account for introntype in cases of ties */
      indel_pos = mismatch_positions_right[righti] + 1;
      left1 = gbuffer[indel_pos];
      left2 = gbuffer[indel_pos+1];
      right2 = gbuffer[indel_pos-indels-2];
      right1 = gbuffer[indel_pos-indels-1];
      introntype = Intron_type(left1,left2,right2,right1,left1,left2,right2,right1,/*cdna_direction*/0);
      intron_level = Intron_level(introntype);
      debug2(printf("  (Case C1) sum %d=%d+%d at indel_pos %d (%c%c-%c%c, type %s).",
		    sum,righti,lefti+1,mismatch_positions_right[righti]+1,
		    left1,left2,right2,right1,Intron_type_string(introntype)));
      if (sum < best_sum ||
	  (sum == best_sum && intron_level > best_intron_level)) {
	if (indel_pos >= min_indel_end_matches && indel_pos <= querylength - min_indel_end_matches) {
	  best_indel_pos = indel_pos;
	  nmismatches_righti = righti;
	  nmismatches_lefti = lefti + 1;
	  debug2(printf("**"));
	  best_sum = sum;
	  *best_introntype = introntype;
	  best_intron_level = intron_level;
	}
      }

    } else {
      debug2(printf("  (Case C1) sum %d=%d+%d at indel_pos %d.",
		    sum,righti,lefti+1,mismatch_positions_right[righti]+1));
      if (sum <= best_sum) {
	indel_pos = mismatch_positions_right[righti] + 1;
	if (indel_pos >= min_indel_end_matches && indel_pos <= querylength - min_indel_end_matches) {
	  best_indel_pos = indel_pos;
	  nmismatches_righti = righti;
	  nmismatches_lefti = lefti + 1;
	  debug2(printf("**"));
	  best_sum = sum;
	}
      }
    }
    righti++;
  }
  debug2(printf("\n"));

  /* Try from other side to see if we missed anything */
  lefti = 0;
  righti = nmismatches_right - 1;

  while (lefti < nmismatches_left) {
    while (righti >= 0 && mismatch_positions_right[righti] < mismatch_positions_left[lefti]) {
      righti--;
    }
    sum = lefti + righti + 1;

    if (-indels >= min_intronlength) {
      /* Account for introntype in cases of ties */
      indel_pos = mismatch_positions_left[lefti];
      left1 = gbuffer[indel_pos];
      left2 = gbuffer[indel_pos+1];
      right2 = gbuffer[indel_pos-indels-2];
      right1 = gbuffer[indel_pos-indels-1];
      introntype = Intron_type(left1,left2,right2,right1,left1,left2,right2,right1,/*cdna_direction*/0);
      intron_level = Intron_level(introntype);
      debug2(printf("  (Case C2) sum %d=%d+%d at indel_pos %d (%c%c-%c%c).",
		    sum,lefti,righti+1,mismatch_positions_left[lefti],
		    gbuffer[indel_pos],gbuffer[indel_pos+1],gbuffer[indel_pos-indels-2],gbuffer[indel_pos-indels-1]));
      if (sum < best_sum ||
	  (sum == best_sum && intron_level > best_intron_level) ||
	  (sum == best_sum && intron_level == best_intron_level && indel_pos < best_indel_pos)) {
	if (indel_pos >= min_indel_end_matches && indel_pos <= querylength - min_indel_end_matches) {
	  best_indel_pos = indel_pos;
	  nmismatches_lefti = lefti;
	  nmismatches_righti = righti + 1;
	  debug2(printf("**"));
	  best_sum = sum;
	  *best_introntype = introntype;
	  best_intron_level = intron_level;
	}
      }

    } else {
      debug2(printf("  (Case C2) sum %d=%d+%d at indel_pos %d.",
		    sum,lefti,righti+1,mismatch_positions_left[lefti]));
      indel_pos = mismatch_positions_left[lefti];
      if (sum < best_sum ||
	  (sum == best_sum && indel_pos < best_indel_pos)) {
	if (indel_pos >= min_indel_end_matches && indel_pos <= querylength - min_indel_end_matches) {
	  best_indel_pos = indel_pos;
	  nmismatches_lefti = lefti;
	  nmismatches_righti = righti + 1;
	  debug2(printf("**"));
	  best_sum = sum;
	}
      }
    }
    lefti++;
  }
  debug2(printf("\n"));

#ifdef DEBUG2
  FREE(gbuffer);
#else
  if (-indels >= min_intronlength) {
    FREE(gbuffer);
  }
#endif

#ifdef HAVE_ALLOCA
  if (querylength <= MAX_STACK_READLENGTH) {
    FREEA(mismatch_positions_left);
    FREEA(mismatch_positions_right);
  } else {
    FREE(mismatch_positions_left);
    FREE(mismatch_positions_right);
  }
#else
  FREE(mismatch_positions_left);
  FREE(mismatch_positions_right);
#endif

  *best_nmismatches_i = nmismatches_lefti;
  *best_nmismatches_j = nmismatches_righti;

  if (best_sum > max_mismatches_allowed) {
    debug2(printf("Returning -1\n"));
    return -1;
#if 0
  } else if (plusp == true) {
    return best_indel_pos;
  } else {
    return querylength - best_indel_pos;
#else
  } else {
    debug2(printf("Returning %d with nmismatches %d+%d\n",
		  best_indel_pos,*best_nmismatches_i,*best_nmismatches_j));
    return best_indel_pos;
#endif
  }
}


/* indels is positive here */
List_T
Indel_solve_middle_insertion (bool *foundp, int *found_score, int *nhits, List_T hits,
			      Univcoord_T left, Chrnum_T chrnum, Univcoord_T chroffset,
			      Univcoord_T chrhigh, Chrpos_T chrlength,
			      int indels, Compress_T query_compress,
			      int querylength, int max_mismatches_allowed,
			      bool plusp, int genestrand, bool sarrayp) {
#ifdef DEBUG2
  int i;
  char *gbuffer;
#endif
  Stage3end_T hit;
  int best_indel_pos, query_indel_pos, indel_pos;
  int nmismatches_left, nmismatches_right;
  int best_sum, sum, nmismatches_lefti, nmismatches_righti, lefti, righti;
  int nmismatches1, nmismatches2;
  int *mismatch_positions_left, *mismatch_positions_right;

#ifdef HAVE_ALLOCA
  if (querylength <= MAX_STACK_READLENGTH) {
    mismatch_positions_left = (int *) ALLOCA(querylength * sizeof(int));
    mismatch_positions_right = (int *) ALLOCA(querylength * sizeof(int));
  } else {
    mismatch_positions_left = (int *) MALLOC(querylength * sizeof(int));
    mismatch_positions_right = (int *) MALLOC(querylength * sizeof(int));
  }
#else
  mismatch_positions_left = (int *) MALLOC(querylength * sizeof(int));
  mismatch_positions_right = (int *) MALLOC(querylength * sizeof(int));
#endif


  *foundp = false;

  /* query has insertion.  Get |indels| less from genome; trim from left. */
  /* left = ptr->diagonal - querylength; */

  assert(indels > 0);
  debug2(gbuffer = (char *) CALLOC(querylength-indels+1,sizeof(char)));
  debug2(Genome_fill_buffer_blocks(left+indels,querylength-indels,gbuffer));
  debug2(printf("solve_middle_indel, plus, insertion: Getting genome at diagonal - querylength %d + indels %d = %llu\n",
		querylength,indels,(unsigned long long) left+indels));
  debug2(printf("g1: %s\n",gbuffer));
  debug2(printf("g2: %s\n",&(gbuffer[indels])));

  /* No need to check chromosome bounds */
  nmismatches_left = Genome_mismatches_left(mismatch_positions_left,max_mismatches_allowed,
					    query_compress,left+indels,/*pos5*/0,/*pos3*/querylength,
					    plusp,genestrand);

  debug2(
	 printf("%d mismatches on left at:",nmismatches_left);
	 for (i = 0; i <= nmismatches_left; i++) {
	   printf(" %d",mismatch_positions_left[i]);
	 }
	 printf("\n");
	 );


  /* No need to check chromosome bounds */
  nmismatches_right = Genome_mismatches_right(mismatch_positions_right,max_mismatches_allowed,
					      query_compress,left,/*pos5*/0,/*pos3*/querylength,
					      plusp,genestrand);

  debug2(
	 printf("%d mismatches on right at:",nmismatches_right);
	 for (i = 0; i <= nmismatches_right; i++) {
	   printf(" %d",mismatch_positions_right[i]);
	 }
	 printf("\n");
	 );

  best_sum = querylength + querylength;

  /* Modeled after end D to get lowest possible coordinate */
  righti = 0;
  lefti = nmismatches_left - 1;
  nmismatches_righti = /*righti*/ 0;
  nmismatches_lefti = /*lefti+1*/ nmismatches_left;

  while (righti < nmismatches_right) {
    while (lefti >= 0 && mismatch_positions_left[lefti] > mismatch_positions_right[righti] - indels) {
      lefti--;
    }
    sum = righti + lefti + 1;
    debug2(printf("  (Case D) sum %d=%d+%d at indel_pos %d.",
		  sum,righti,lefti+1,mismatch_positions_right[righti]-indels+1));
    if (sum <= best_sum) {
      indel_pos = mismatch_positions_right[righti] - indels + 1;
      if (indel_pos >= min_indel_end_matches && indel_pos + indels <= querylength - min_indel_end_matches &&
	  /* Require more matches on both ends than insertion length */ indel_pos >= indels && indel_pos + indels <= querylength - indels) {
	best_indel_pos = indel_pos;
	nmismatches_righti = righti;
	nmismatches_lefti = lefti + 1;
	debug2(printf("**"));
	best_sum = sum;
      }
    }
    righti++;
  }
  debug2(printf("\n"));


  /* Try from other side to see if we missed anything */
  lefti = 0;
  righti = nmismatches_right - 1;

  while (lefti < nmismatches_left) {
    while (righti >= 0 && mismatch_positions_right[righti] < mismatch_positions_left[lefti] + indels) {
      righti--;
    }
    sum = lefti + righti + 1;
    debug2(printf("  (Case D2) sum %d=%d+%d at indel_pos %d.",
		  sum,lefti,righti+1,mismatch_positions_left[lefti]));
    if (sum < best_sum) {
      indel_pos = mismatch_positions_left[lefti];
      if (indel_pos >= min_indel_end_matches && indel_pos + indels <= querylength - min_indel_end_matches &&
	  /* Require more matches on both ends than insertion length */ indel_pos >= indels && indel_pos + indels <= querylength - indels) {
	best_indel_pos = indel_pos;
	nmismatches_righti = righti + 1;
	nmismatches_lefti = lefti;
	debug2(printf("**"));
	best_sum = sum;
      }
    } else if (sum == best_sum) {
      indel_pos = mismatch_positions_left[lefti];
      if (indel_pos < best_indel_pos) {
	if (indel_pos >= min_indel_end_matches && indel_pos + indels <= querylength - min_indel_end_matches &&
	  /* Require more matches on both ends than insertion length */ indel_pos >= indels && indel_pos + indels <= querylength - indels) {
	  best_indel_pos = indel_pos;
	  nmismatches_righti = righti + 1;
	  nmismatches_lefti = lefti;
	  debug2(printf("**"));
	  /* best_sum = sum; */
	}
      }
    }
    lefti++;
  }
  debug2(printf("\n"));

#ifdef HAVE_ALLOCA
  if (querylength <= MAX_STACK_READLENGTH) {
    FREEA(mismatch_positions_left);
    FREEA(mismatch_positions_right);
  } else {
    FREE(mismatch_positions_left);
    FREE(mismatch_positions_right);
  }
#else
  FREE(mismatch_positions_left);
  FREE(mismatch_positions_right);
#endif

  if (best_sum <= max_mismatches_allowed) {
    if (plusp == true) {
      query_indel_pos = best_indel_pos;
      nmismatches1 = nmismatches_lefti;
      nmismatches2 = nmismatches_righti;
    } else {
      query_indel_pos = querylength - best_indel_pos - indels;
      nmismatches1 = nmismatches_righti;
      nmismatches2 = nmismatches_lefti;
    }

    /* genomiclength = querylength-indels; */
    if ((hit = Stage3end_new_insertion(&(*found_score),indels,query_indel_pos,
				       nmismatches1,nmismatches2,/*left*/left+indels,
				       query_compress,querylength,plusp,genestrand,
				       chrnum,chroffset,chrhigh,chrlength,
				       indel_penalty_middle,sarrayp)) != NULL) {
      debug2(printf("successful insertion with %d=%d+%d mismatches and indel_pos at %d\n",
		    sum,nmismatches_lefti,nmismatches_righti,best_indel_pos));
      /* ptr->usedp = ptr2->usedp = true; */
      *foundp = true;
      *nhits += 1;
      hits = List_push(hits,(void *) hit);
    }
  }

  return hits;
}



/* indels is negative here */
List_T
Indel_solve_middle_deletion (bool *foundp, int *found_score, int *nhits, List_T hits,
			     Univcoord_T left, Chrnum_T chrnum, Univcoord_T chroffset,
			     Univcoord_T chrhigh, Chrpos_T chrlength,
			     int indels, Compress_T query_compress, int querylength,
			     int max_mismatches_allowed,
			     bool plusp, int genestrand, bool sarrayp) {
#ifdef DEBUG2
  int i;
  char *gbuffer;
#endif
  Stage3end_T hit;
  int best_indel_pos, query_indel_pos, indel_pos;
  int nmismatches_left, nmismatches_right;
  int best_sum, sum, nmismatches_lefti, nmismatches_righti, lefti, righti;
  int nmismatches1, nmismatches2;
  int *mismatch_positions_left, *mismatch_positions_right;

#ifdef HAVE_ALLOCA
  if (querylength <= MAX_STACK_READLENGTH) {
    mismatch_positions_left = (int *) ALLOCA(querylength * sizeof(int));
    mismatch_positions_right = (int *) ALLOCA(querylength * sizeof(int));
  } else {
    mismatch_positions_left = (int *) MALLOC(querylength * sizeof(int));
    mismatch_positions_right = (int *) MALLOC(querylength * sizeof(int));
  }
#else
  mismatch_positions_left = (int *) MALLOC(querylength * sizeof(int));
  mismatch_positions_right = (int *) MALLOC(querylength * sizeof(int));
#endif


  *foundp = false;

  /* query has deletion.  Get |indels| more from genome; add to right. */
  /* left = ptr->diagonal - querylength; */

  assert(indels < 0);
  debug2(gbuffer = (char *) CALLOC(querylength-indels+1,sizeof(char)));
  debug2(Genome_fill_buffer_blocks(left,querylength-indels,gbuffer));
  debug2(printf("solve_middle_indel, plus, deletion (indels %d): Getting genome at diagonal - querylength %d = %llu\n",
		indels,querylength,(unsigned long long) left));
  debug2(printf("g1: %s\n",gbuffer));
  debug2(printf("g2: %s\n",&(gbuffer[-indels])));
  debug2(FREE(gbuffer));

  /* No need to check chromosome bounds */
  nmismatches_left = Genome_mismatches_left(mismatch_positions_left,max_mismatches_allowed,
					    query_compress,left,/*pos5*/0,/*pos3*/querylength,
					    plusp,genestrand);

  debug2(
	 printf("%d mismatches on left at:",nmismatches_left);
	 for (i = 0; i <= nmismatches_left; i++) {
	   printf(" %d",mismatch_positions_left[i]);
	 }
	 printf("\n");
	 );

  /* No need to check chromosome bounds */
  nmismatches_right = Genome_mismatches_right(mismatch_positions_right,max_mismatches_allowed,
					      query_compress,left-indels,/*pos5*/0,/*pos3*/querylength,
					      plusp,genestrand);

  debug2(
	 printf("%d mismatches on right at:",nmismatches_right);
	 for (i = 0; i <= nmismatches_right; i++) {
	   printf(" %d",mismatch_positions_right[i]);
	 }
	 printf("\n");
	 );

  best_sum = querylength + querylength;

  /* Modeled after end C to get lowest possible coordinate */
  righti = 0;
  lefti = nmismatches_left - 1;
  nmismatches_righti = /*righti*/ 0;
  nmismatches_lefti = /*lefti+1*/ nmismatches_left;

  while (righti < nmismatches_right) {
    while (lefti >= 0 && mismatch_positions_left[lefti] > mismatch_positions_right[righti]) {
      lefti--;
    }
    sum = righti + lefti + 1;
    debug2(printf("  (Case C1) sum %d=%d+%d at indel_pos %d.",
		  sum,righti,lefti+1,mismatch_positions_right[righti]+1));
    if (sum <= best_sum) {
      indel_pos = mismatch_positions_right[righti] + 1;
      if (indel_pos >= min_indel_end_matches && indel_pos <= querylength - min_indel_end_matches) {
	best_indel_pos = indel_pos;
	nmismatches_righti = righti;
	nmismatches_lefti = lefti + 1;
	debug2(printf("**"));
	best_sum = sum;
      }
    }
    righti++;
  }
  debug2(printf("\n"));

  /* Try from other side to see if we missed anything */
  lefti = 0;
  righti = nmismatches_right - 1;

  while (lefti < nmismatches_left) {
    while (righti >= 0 && mismatch_positions_right[righti] < mismatch_positions_left[lefti]) {
      righti--;
    }
    sum = lefti + righti + 1;
    debug2(printf("  (Case C2) sum %d=%d+%d at indel_pos %d.",
		  sum,lefti,righti+1,mismatch_positions_left[lefti]));
    if (sum < best_sum) {
      indel_pos = mismatch_positions_left[lefti];
      if (indel_pos >= min_indel_end_matches && indel_pos <= querylength - min_indel_end_matches) {
	best_indel_pos = indel_pos;
	nmismatches_lefti = lefti;
	nmismatches_righti = righti + 1;
	debug2(printf("**"));
	best_sum = sum;
      }
    } else if (sum == best_sum) {
      indel_pos = mismatch_positions_left[lefti];
      if (indel_pos < best_indel_pos) {
	if (indel_pos >= min_indel_end_matches && indel_pos <= querylength - min_indel_end_matches) {
	  best_indel_pos = indel_pos;
	  nmismatches_lefti = lefti;
	  nmismatches_righti = righti + 1;
	  debug2(printf("**"));
	  /* best_sum = sum; */
	}
      }
    }
    lefti++;
  }
  debug2(printf("\n"));

#ifdef HAVE_ALLOCA
  if (querylength <= MAX_STACK_READLENGTH) {
    FREEA(mismatch_positions_left);
    FREEA(mismatch_positions_right);
  } else {
    FREE(mismatch_positions_left);
    FREE(mismatch_positions_right);
  }
#else
  FREE(mismatch_positions_left);
  FREE(mismatch_positions_right);
#endif

  if (best_sum <= max_mismatches_allowed) {
    if (plusp == true) {
      query_indel_pos = best_indel_pos;
      nmismatches1 = nmismatches_lefti;
      nmismatches2 = nmismatches_righti;
    } else {
      query_indel_pos = querylength - best_indel_pos;
      nmismatches1 = nmismatches_righti;
      nmismatches2 = nmismatches_lefti;
    }

    /* genomiclength = querylength-indels; */
    if ((hit = Stage3end_new_deletion(&(*found_score),-indels,query_indel_pos,
				      nmismatches1,nmismatches2,left,
				      query_compress,querylength,plusp,genestrand,
				      chrnum,chroffset,chrhigh,chrlength,
				      indel_penalty_middle,sarrayp)) != NULL) {
      debug2(printf("successful middle deletion with %d=%d+%d mismatches and indel_pos at %d and nindels %d\n",
		    best_sum,nmismatches_lefti,nmismatches_righti,best_indel_pos,-indels));
      *foundp = true;
      /* ptr->usedp = ptr2->usedp = true; */
      *nhits += 1;
      hits = List_push(hits,(void *) hit);
    }
  }

  return hits;
}