@@ -1,4 +1,4 @@

-static char rcsid[] = "$Id: bamtally.c 135743 2014-05-09 16:58:57Z twu $";

+static char rcsid[] = "$Id: bamtally.c 136268 2014-05-14 22:04:03Z michafla $";

 #ifdef HAVE_CONFIG_H

 #include <config.h>

 #endif

@@ -5037,7 +5037,7 @@ Bamtally_iit (Bamreader_T bamreader, char *desired_chr, char *bam_lacks_chr,

 		   /*quality_score_adj*/0,min_depth,variant_strands,

 		   /*genomic_diff_p*/false,/*signed_counts_p*/false,ignore_query_Ns_p,

 		   print_indels_p,/*print_totals_p*/false,/*print_cycles_p*/false,

-		   /*print_quality_scores_p*/false,/*print_mapq_scores_p*/false,/*print_xs_scores_p*/print_xs_scores_p,

+		   /*print_quality_scores_p*/false,/*print_mapq_scores_p*/false,print_xs_scores_p,

 		   /*want_genotypes_p*/false,verbosep,readlevel_p,max_softclip,

 		   print_noncovered_p,/*bamfile*/NULL);

       /* Reverse lists so we can specify presortedp == true */

@@ -5087,7 +5087,7 @@ Bamtally_iit (Bamreader_T bamreader, char *desired_chr, char *bam_lacks_chr,

                    /*genomic_diff_p*/false,/*signed_counts_p*/false,

                    ignore_query_Ns_p, print_indels_p,/*print_totals_p*/false,

                    /*print_cycles_p*/false,

-                   /*print_quality_scores_p*/false,/*print_mapq_scores_p*/false,/*print_xs_scores_p*/print_xs_scores_p,

+                   /*print_quality_scores_p*/false,/*print_mapq_scores_p*/false,print_xs_scores_p,

                    /*want_genotypes_p*/false,verbosep,readlevel_p,max_softclip,

 		   print_noncovered_p,/*bamfile*/NULL);

     }

@@ -1,4 +1,4 @@

-static char rcsid[] = "$Id: dynprog.c 133253 2014-04-15 18:56:11Z twu $";

+static char rcsid[] = "$Id: dynprog.c 136514 2014-05-16 17:59:19Z twu $";

 #ifdef HAVE_CONFIG_H

 #include <config.h>

 #endif

@@ -75,11 +75,7 @@ static char rcsid[] = "$Id: dynprog.c 133253 2014-04-15 18:56:11Z twu $";

 #define FULLMATCH 3

 #define HALFMATCH 1

-#ifdef PMAP

 #define AMBIGUOUS 0

-#else

-#define AMBIGUOUS -1

-#endif

 /* These values were set to -5, -4, -3, but this led to chopped ends

@@ -764,6 +760,7 @@ Pairdistance_T **pairdistance_array[NMISMATCHTYPES];

 Pairdistance_T **pairdistance_array_plus_128[NMISMATCHTYPES];

 #endif

 bool **consistent_array;

+int *nt_to_int_array;

 bool

@@ -854,6 +851,16 @@ Dynprog_init (Mode_T mode) {

   int i, j, ptr;

   int c, c1, c2;

+  nt_to_int_array = (int *) CALLOC(128,sizeof(int));

+  for (j = 0; j < 128; j++) {

+    nt_to_int_array[j] = 4;

+  }

+  nt_to_int_array['A'] = nt_to_int_array['a'] = 0;

+  nt_to_int_array['C'] = nt_to_int_array['c'] = 1;

+  nt_to_int_array['G'] = nt_to_int_array['g'] = 2;

+  nt_to_int_array['T'] = nt_to_int_array['t'] = 3;

+

+

   consistent_array = (bool **) CALLOC(128,sizeof(bool *));

   consistent_array[0] = (bool *) CALLOC(128*128,sizeof(bool));

   ptr = 0;

@@ -898,6 +905,11 @@ Dynprog_init (Mode_T mode) {

   }

 #endif

+  for (c = 'A'; c < 'Z'; c++) {

+    permute_cases(c,c,FULLMATCH);

+  }

+

+  /* Exceptions */

   permute_cases('U','T',FULLMATCH);

   permute_cases('R','A',HALFMATCH);

@@ -944,16 +956,16 @@ Dynprog_init (Mode_T mode) {

   permute_cases('X','A',AMBIGUOUS);

   permute_cases('X','G',AMBIGUOUS);

+  permute_cases('N','N',AMBIGUOUS); /* Needed to start dynprog procedures with 0 at (0,0) */

+  permute_cases('X','X',AMBIGUOUS); /* Needed to start dynprog procedures with 0 at (0,0) */

+

+

   if (mode == CMET_STRANDED || mode == CMET_NONSTRANDED) {

     /* Query-T can match Genomic-C */

     permute_cases_oneway('T','C',FULLMATCH);

     permute_cases_oneway('A','G',FULLMATCH);

   }

-  for (c = 'A'; c < 'Z'; c++) {

-    permute_cases(c,c,FULLMATCH);

-  }

-

 #ifndef HAVE_SSE4_1

 #ifdef PREUC

   for (i = 0; i < NMISMATCHTYPES; i++) {

@@ -1010,6 +1022,7 @@ Dynprog_term (void) {

   */

   FREE(consistent_array[0]);

   FREE(consistent_array);

+  FREE(nt_to_int_array);

   return;

 }

@@ -1678,274 +1691,3 @@ Dynprog_traceback_std (List_T pairs, int *nmatches, int *nmismatches, int *nopen

 #endif

-char *

-Dynprog_tokens_string (List_T tokens) {

-  char *string, *token, *q;

-  int stringlength = 0;

-  List_T p;

-

-  for (p = tokens; p != NULL; p = List_next(p)) {

-    token = (char *) List_head(p);

-    stringlength += strlen(token);

-  }

-  string = (char *) MALLOC((stringlength+1) * sizeof(char));

-

-  q = string;

-  for (p = tokens; p != NULL; p = List_next(p)) {

-    token = (char *) List_head(p);

-    strcpy(q,token);

-    q += strlen(token);

-  }

-  *q = '\0';

-  return string;

-}

-

-

-void

-Dynprog_tokens_free (List_T *tokens) {

-  List_T p;

-  char *token;

-

-  for (p = *tokens; p != NULL; p = List_next(p)) {

-    token = (char *) List_head(p);

-    FREE(token);

-  }

-  List_free(&(*tokens));

-

-  return;

-}

-

-static List_T

-push_token (bool *startp, List_T tokens, char *token) {

-  char *copy;

-

-  copy = (char *) CALLOC(strlen(token)+1,sizeof(char));

-  strcpy(copy,token);

-

-  *startp = false;

-  return List_push(tokens,(void *) copy);

-}

-

-

-

-char *

-Dynprog_cigar_std (int *finalc, Direction32_T **directions_nogap, Direction32_T **directions_Egap, Direction32_T **directions_Fgap,

-		   int r, int c, char *rsequence, char *gsequence, char *gsequence_alt,

-		   char *nindels, int queryoffset, int genomeoffset, bool revp,

-		   Univcoord_T chroffset, Univcoord_T chrhigh) {

-  char *cigar;

-  List_T tokens = NULL;

-  char token[10];

-  int Mlength = 0, Ilength = 0, Mlength_postins = 0;

-  int insertion_width = 0;

-  bool startp = true;

-

-  int dist;

-  Direction32_T dir;

-

-  debug(printf("Starting Dynprog_cigar_std at r=%d,c=%d (roffset=%d, goffset=%d)\n",r,c,queryoffset,genomeoffset));

-

-#if 0

-  /* Handle initial indel */

-  if ((dir = directions_nogap[c][r]) == DIAG) {

-    /* Not an indel.  Do nothing. */

-

-  } else if (dir == HORIZ) {

-    dist = 1;

-    while (c > 1 && directions_Egap[c][r] != DIAG) {

-      dist++;

-      c--;

-    }

-    c--;

-    /* dir = directions_nogap[c][r]; */

-

-    sprintf(token,"0M");

-    tokens = push_token(&startp,tokens,token);

-    sprintf(token,"%dD",dist);

-    tokens = push_token(&startp,tokens,token);

-

-  } else {

-    /* Must be VERT */

-    dist = 1;

-    while (r > 1 && directions_Fgap[c][r] != DIAG) {

-      dist++;

-      r--;

-    }

-    r--;

-    /* dir = directions_nogap[c][r]; */

-

-    sprintf(token,"0M");

-    tokens = push_token(&startp,tokens,token);

-    sprintf(token,"%dI",dist);

-    tokens = push_token(&startp,tokens,token);

-  }

-#endif

-

-  *finalc = c;

-  while (c > 0) {  /* dir != STOP */

-    if (r == 0) {

-      /* Ignore gap in row 0 */

-      if (Ilength > 0 && Ilength != insertion_width) {

-	if (Mlength_postins > 0) {

-	  sprintf(token,"%dM",Mlength_postins);

-	  tokens = push_token(&startp,tokens,token);

-	}

-	sprintf(token,"%dS",Ilength);

-	tokens = push_token(&startp,tokens,token);

-	if (Mlength > 0) {

-	  sprintf(token,"%dM",Mlength);

-	  tokens = push_token(&startp,tokens,token);

-	}

-      } else {

-	if (Ilength > 0) {

-	  sprintf(token,"%dM",Mlength_postins);

-	  tokens = push_token(&startp,tokens,token);

-	  sprintf(token,"%dI",Ilength);

-	  tokens = push_token(&startp,tokens,token);

-	}

-	if (1 || Mlength > 0) {

-	  sprintf(token,"%dM",Mlength);

-	  tokens = push_token(&startp,tokens,token);

-	}

-      }

-

-      *finalc = c;

-      tokens = List_reverse(tokens);

-      cigar = Dynprog_tokens_string(List_reverse(tokens));

-      Dynprog_tokens_free(&tokens);

-

-      return cigar;

-

-    } else if ((dir = directions_nogap[c][r]) == DIAG) {

-      /* printf("At r %d, c %d (%c), dir is DIAG, nindels %d\n",r,c,gsequence[c],nindels[c-1]); */

-

-      if (nindels[c-1] < 0) {

-	/* Genome modified with a deletion */

-	if (startp == true && Ilength > 0 && Ilength != insertion_width) {

-	  /* Incomplete insertion at beginning.  Mlength_postins should be 0. */

-	  sprintf(token,"%dS",Ilength);

-	  tokens = push_token(&startp,tokens,token);

-	  sprintf(token,"%dM",Mlength);

-	  tokens = push_token(&startp,tokens,token);

-	} else {

-	  if (Ilength > 0) {

-	    sprintf(token,"%dM",Mlength_postins);

-	    tokens = push_token(&startp,tokens,token);

-	    sprintf(token,"%dI",Ilength);

-	    tokens = push_token(&startp,tokens,token);

-	  }

-	  if (1 || Mlength > 0) {

-	    sprintf(token,"%dM",Mlength);

-	    tokens = push_token(&startp,tokens,token);

-	  }

-	}

-	Mlength = 1;

-

-	sprintf(token,"%dD",-nindels[c-1]);

-	tokens = push_token(&startp,tokens,token);

-

-      } else if (nindels[c-1] > 0) {

-	insertion_width = nindels[c-1];

-	if (Ilength == 0) {

-	  Mlength_postins = Mlength; /* Cannot push M yet, since I could change to S */

-	  Mlength = 0;

-	}

-	Ilength += 1;

-	/* printf("Incrementing Ilength to be %d\n",Ilength); */

-

-      } else {

-	Mlength += 1;

-	/* printf("Incrementing Mlength to be %d\n",Mlength); */

-      }

-      r--; c--;

-

-    } else if (dir == HORIZ) {

-      /* printf("At r %d, c %d, dir is HORIZ\n",r,c); */

-

-      if (startp == true && Ilength > 0 && Ilength != insertion_width) {

-	/* Incomplete insertion at beginning.  Mlength_postins should be 0. */

-	sprintf(token,"%dS",Ilength);

-	tokens = push_token(&startp,tokens,token);

-	sprintf(token,"%dM",Mlength);

-	tokens = push_token(&startp,tokens,token);

-      } else {

-	if (Ilength > 0) {

-	  sprintf(token,"%dM",Mlength_postins);

-	  tokens = push_token(&startp,tokens,token);

-	  sprintf(token,"%dI",Ilength);

-	  tokens = push_token(&startp,tokens,token);

-	}

-	if (1 || Mlength > 0) {

-	  sprintf(token,"%dM",Mlength);

-	  tokens = push_token(&startp,tokens,token);

-	}

-      }

-      Mlength = 0;

-

-      dist = 1;

-      while (c > 1 && directions_Egap[c][r] != DIAG) {

-	dist++;

-	c--;

-      }

-      c--;

-      /* dir = directions_nogap[c][r]; */

-

-      sprintf(token,"%dD",dist);

-      tokens = push_token(&startp,tokens,token);

-

-    } else {

-      /* Must be VERT */

-      /* printf("At r %d, c %d, dir is VERT\n",r,c); */

-      if (Ilength == 0) {

-	if (1 || Mlength > 0) {

-	  sprintf(token,"%dM",Mlength);

-	  tokens = push_token(&startp,tokens,token);

-	  Mlength = 0;

-	}

-      }

-

-      Ilength += 1;

-      while (r > 1 && directions_Fgap[c][r] != DIAG) {

-	Ilength += 1;

-	r--;

-      }

-      r--;

-      /* dir = directions_nogap[c][r]; */

-

-      sprintf(token,"%dI",Ilength);

-      tokens = push_token(&startp,tokens,token);

-      Ilength = 0;

-    }

-  }

-

-  if (Ilength > 0 && Ilength != insertion_width) {

-    if (Mlength_postins > 0) {

-      sprintf(token,"%dM",Mlength_postins);

-      tokens = push_token(&startp,tokens,token);

-    }

-    sprintf(token,"%dS",Ilength);

-    tokens = push_token(&startp,tokens,token);

-    if (Mlength > 0) {

-      sprintf(token,"%dM",Mlength);

-      tokens = push_token(&startp,tokens,token);

-    }

-  } else {

-    if (Ilength > 0) {

-      sprintf(token,"%dM",Mlength_postins);

-      tokens = push_token(&startp,tokens,token);

-      sprintf(token,"%dI",Ilength);

-      tokens = push_token(&startp,tokens,token);

-    }

-    if (1 || Mlength > 0) {

-      sprintf(token,"%dM",Mlength);

-      tokens = push_token(&startp,tokens,token);

-    }

-  }

-

-  *finalc = c;

-  tokens = List_reverse(tokens);

-  cigar = Dynprog_tokens_string(List_reverse(tokens));

-  Dynprog_tokens_free(&tokens);

-

-  return cigar;

-}

@@ -1,4 +1,4 @@

-/* $Id: dynprog.h 133170 2014-04-14 23:43:26Z twu $ */

+/* $Id: dynprog.h 136514 2014-05-16 17:59:19Z twu $ */

 #ifndef DYNPROG_INCLUDED

 #define DYNPROG_INCLUDED

@@ -54,9 +54,11 @@ typedef int Direction32_T;

 #define DIAG 0			/* Pre-dominant case.  Directions_alloc clears to this value. */

 #define NEG_INFINITY_8 (-128)

+#define POS_INFINITY_8 (127)

 #define MAX_CHAR (127)

 #define NEG_INFINITY_16 (-32768)

+#define POS_INFINITY_16 (32767)

 #define MAX_SHORT (32767)

 /* We can allow -128 and -32768 for NEG_INFINITY in SIMD procedures,

@@ -85,6 +87,7 @@ extern Pairdistance_T **pairdistance_array[NMISMATCHTYPES];

 extern Pairdistance_T **pairdistance_array_plus_128[NMISMATCHTYPES];

 #endif

 extern bool **consistent_array;

+extern int *nt_to_int_array;

 struct Space_single_T {

@@ -171,14 +174,6 @@ Dynprog_tokens_string (List_T tokens);

 extern void

 Dynprog_tokens_free (List_T *tokens);

-extern char *

-Dynprog_cigar_std (int *finalc, Direction32_T **directions_nogap, Direction32_T **directions_Egap,

-		   Direction32_T **directions_Fgap,

-		   int r, int c, char *rsequence, char *gsequence, char *gsequence_alt,

-		   char *nindels, int queryoffset, int genomeoffset, bool revp,

-		   Univcoord_T chroffset, Univcoord_T chrhigh);

-

-

 #undef T

 #endif

@@ -1,4 +1,4 @@

-static char rcsid[] = "$Id: dynprog_simd.c 133254 2014-04-15 18:56:28Z twu $";

+static char rcsid[] = "$Id: dynprog_simd.c 136516 2014-05-16 18:00:16Z twu $";

 #ifdef HAVE_CONFIG_H

 #include <config.h>

 #endif

@@ -17,6 +17,18 @@ static char rcsid[] = "$Id: dynprog_simd.c 133254 2014-04-15 18:56:28Z twu $";

 #endif

 #include "mem.h"

+#include "assert.h"

+

+

+#define LAZY_INDEL 1		/* Don't advance to next coordinate on final indel, since could go over chromosome bounds. */

+

+/* Row 0 and column 0 directions */

+/* Was useful in finding a saturation bug, but can fail because of saturation */

+#ifdef CHECK1

+#define check1(x) x

+#else

+#define check1(x)

+#endif

 #ifdef DEBUG

@@ -31,13 +43,19 @@ static char rcsid[] = "$Id: dynprog_simd.c 133254 2014-04-15 18:56:28Z twu $";

 #define debug2(x)

 #endif

-/* F loop */

+/* Fgap */

 #ifdef DEBUG3

 #define debug3(x) x

 #else

 #define debug3(x)

 #endif

+#ifdef DEBUG8

+#define debug8(x) x

+#else

+#define debug8(x)

+#endif

+

 #ifdef DEBUG14

 #define debug14(x) x

 #else

@@ -52,13 +70,13 @@ static char rcsid[] = "$Id: dynprog_simd.c 133254 2014-04-15 18:56:28Z twu $";

-/************************************************************************

- *   Debugging procedures

- ************************************************************************/

-

+#include "complement.h"

 #define NEG_INFINITY_DISPLAY -99

+/************************************************************************

+ *   Debugging procedures

+ ************************************************************************/

 #ifdef DEBUG15

 /* For debugging of SIMD procedures*/

@@ -87,20 +105,27 @@ print_vector_16 (__m128i x, int r, int c, char *label) {

 #endif

-#if defined(DEBUG2) || defined(DEBUG14)

+#if defined(DEBUG14) || defined(DEBUG2)

 static void

 Matrix8_print (Score8_T **matrix, int rlength, int glength, char *rsequence,

 	       char *gsequence, char *gsequencealt,

-#ifdef DEBUG14

-	       int goffset, Univcoord_T chroffset, Univcoord_T chrhigh, bool watsonp,

-#endif

 	       bool revp, int lband, int uband) {

   int i, j;

-  char g_alt;

+  char g, g_alt;

   _mm_lfence();

+  /* j */

+  printf("   ");		/* For i */

+  printf("  ");

+  for (j = 0; j <= glength; ++j) {

+    printf(" %2d ",j);

+  }

+  printf("\n");

+

+

   if (gsequence) {

+    printf("   ");		/* For i */

     printf("  ");

     for (j = 0; j <= glength; ++j) {

       if (j == 0) {

@@ -112,23 +137,27 @@ Matrix8_print (Score8_T **matrix, int rlength, int glength, char *rsequence,

     printf("\n");

   }

-#if 0

-  printf("  ");

-  for (j = 0; j <= glength; ++j) {

-    if (j == 0) {

-      printf("    ");

-    } else {

-      if (revp == false) {

-	printf("  %c ",get_genomic_nt(&g_alt,goffset+j-1,chroffset,chrhigh,watsonp));

+  if (gsequencealt != gsequence) {

+    printf("   ");		/* For i */

+    printf("  ");

+    for (j = 0; j <= glength; ++j) {

+      if (j == 0) {

+	printf("    ");

       } else {

-	printf("  %c ",get_genomic_nt(&g_alt,goffset+1-j,chroffset,chrhigh,watsonp));

+	g = revp ? gsequence[-j+1] : gsequence[j-1];

+	g_alt = revp ? gsequencealt[-j+1] : gsequencealt[j-1];

+	if (g == g_alt) {

+	  printf("  %c ",' ');

+	} else {

+	  printf("  %c ",g_alt);

+	}

       }

     }

+    printf("\n");

   }

-  printf("\n");

-#endif

   for (i = 0; i <= rlength; ++i) {

+    printf("%2d ",i);

     if (i == 0) {

       printf("  ");

     } else {

@@ -152,19 +181,114 @@ Matrix8_print (Score8_T **matrix, int rlength, int glength, char *rsequence,

   return;

 }

+static void

+Matrix8_print_ud (Score8_T **matrix, int rlength, int glength, char *rsequence,

+		  char *gsequence, char *gsequencealt,

+		  bool revp, int band, bool upperp) {

+  int i, j;

+  char g, g_alt;

+

+  _mm_lfence();

+

+  /* j */

+  printf("   ");		/* For i */

+  printf("  ");

+  for (j = 0; j <= glength; ++j) {

+    printf(" %2d ",j);

+  }

+  printf("\n");

+

+  if (gsequence) {

+    printf("   ");		/* For i */

+    printf("  ");

+    for (j = 0; j <= glength; ++j) {

+      if (j == 0) {

+	printf("    ");

+      } else {

+	printf("  %c ",revp ? gsequence[-j+1] : gsequence[j-1]);

+      }

+    }

+    printf("\n");

+  }

+

+  if (gsequencealt != gsequence) {

+    printf("   ");		/* For i */

+    printf("  ");

+    for (j = 0; j <= glength; ++j) {

+      if (j == 0) {

+	printf("    ");

+      } else {

+	g = revp ? gsequence[-j+1] : gsequence[j-1];

+	g_alt = revp ? gsequencealt[-j+1] : gsequencealt[j-1];

+	if (g == g_alt) {

+	  printf("  %c ",' ');

+	} else {

+	  printf("  %c ",g_alt);

+	}

+      }

+    }

+    printf("\n");

+  }

+

+  for (i = 0; i <= rlength; ++i) {

+    printf("%2d ",i);

+    if (i == 0) {

+      printf("  ");

+    } else {

+      printf("%c ",revp ? rsequence[-i+1] : rsequence[i-1]);

+    }

+    if (upperp == true) {

+      for (j = 0; j <= glength; ++j) {

+	if (j < i) {

+	  printf("  . ");

+	} else if (j > i + band) {

+	  printf("  . ");

+	} else if (matrix[j][i] < NEG_INFINITY_DISPLAY) {

+	  printf("%3d ",NEG_INFINITY_DISPLAY);

+	} else {

+	  printf("%3d ",matrix[j][i]);

+	}

+      }

+    } else {

+      for (j = 0; j <= glength; ++j) {

+	if (i < j) {

+	  printf("  . ");

+	} else if (i > j + band) {

+	  printf("  . ");

+	} else if (matrix[i][j] < NEG_INFINITY_DISPLAY) {

+	  printf("%3d ",NEG_INFINITY_DISPLAY);

+	} else {

+	  printf("%3d ",matrix[i][j]);

+	}

+      }

+    }

+    printf("\n");

+  }

+  printf("\n");

+

+  return;

+}

+

+

 static void

 Matrix16_print (Score16_T **matrix, int rlength, int glength, char *rsequence,

 		char *gsequence, char *gsequencealt,

-#ifdef DEBUG14

-		int goffset, Univcoord_T chroffset, Univcoord_T chrhigh, bool watsonp,

-#endif

 		bool revp, int lband, int uband) {

   int i, j;

-  char g_alt;

+  char g, g_alt;

   _mm_lfence();

+  /* j */

+  printf("   ");		/* For i */

+  printf("  ");

+  for (j = 0; j <= glength; ++j) {

+    printf(" %2d ",j);

+  }

+  printf("\n");

+

   if (gsequence) {

+    printf("   ");		/* For i */

     printf("  ");

     for (j = 0; j <= glength; ++j) {

       if (j == 0) {

@@ -176,23 +300,27 @@ Matrix16_print (Score16_T **matrix, int rlength, int glength, char *rsequence,

     printf("\n");

   }

-#if 0

-  printf("  ");

-  for (j = 0; j <= glength; ++j) {

-    if (j == 0) {

-      printf("    ");

-    } else {

-      if (revp == false) {

-	printf("  %c ",get_genomic_nt(&g_alt,goffset+j-1,chroffset,chrhigh,watsonp));

+  if (gsequencealt != gsequence) {

+    printf("   ");		/* For i */

+    printf("  ");

+    for (j = 0; j <= glength; ++j) {

+      if (j == 0) {

+	printf("    ");

       } else {

-	printf("  %c ",get_genomic_nt(&g_alt,goffset+1-j,chroffset,chrhigh,watsonp));

+	g = revp ? gsequence[-j+1] : gsequence[j-1];

+	g_alt = revp ? gsequencealt[-j+1] : gsequencealt[j-1];

+	if (g == g_alt) {

+	  printf("  %c ",' ');

+	} else {

+	  printf("  %c ",g_alt);

+	}

       }

     }

+    printf("\n");

   }

-  printf("\n");

-#endif

   for (i = 0; i <= rlength; ++i) {

+    printf("%2d ",i);

     if (i == 0) {

       printf("  ");

     } else {

@@ -216,105 +344,115 @@ Matrix16_print (Score16_T **matrix, int rlength, int glength, char *rsequence,

   return;

 }

-#endif

-

-#if defined(DEBUG2) || defined(DEBUG14)

 static void

-Directions8_print (Direction8_T **directions_nogap, Direction8_T **directions_Egap, Direction8_T **directions_Fgap,

-		   int rlength, int glength, char *rsequence, char *gsequence, char *gsequence_alt,

-#ifdef DEBUG14

-		   int goffset, Univcoord_T chroffset, Univcoord_T chrhigh, bool watsonp,

-#endif

-		   bool revp, int lband, int uband) {

+Matrix16_print_ud (Score16_T **matrix, int rlength, int glength, char *rsequence,

+		   char *gsequence, char *gsequencealt,

+		   bool revp, int band, bool upperp) {

   int i, j;

-  char g_alt;

+  char g, g_alt;

   _mm_lfence();

+  /* j */

+  printf("   ");		/* For i */

+  printf("  ");

+  for (j = 0; j <= glength; ++j) {

+    printf(" %2d ",j);

+  }

+  printf("\n");

+

   if (gsequence) {

+    printf("   ");		/* For i */

     printf("  ");

     for (j = 0; j <= glength; ++j) {

       if (j == 0) {

-	printf("      ");

+	printf("    ");

       } else {

-	printf("  %c   ",revp ? gsequence[-j+1] : gsequence[j-1]);

+	printf("  %c ",revp ? gsequence[-j+1] : gsequence[j-1]);

       }

     }

     printf("\n");

   }

-#if 0

-  printf("  ");

-  for (j = 0; j <= glength; ++j) {

-    if (j == 0) {

-      printf("      ");

-    } else {

-      if (revp == false) {

-	printf("  %c   ",get_genomic_nt(&g_alt,goffset+j-1,chroffset,chrhigh,watsonp));

+  if (gsequencealt != gsequence) {

+    printf("   ");		/* For i */

+    printf("  ");

+    for (j = 0; j <= glength; ++j) {

+      if (j == 0) {

+	printf("    ");

       } else {

-	printf("  %c   ",get_genomic_nt(&g_alt,goffset+1-j,chroffset,chrhigh,watsonp));

+	g = revp ? gsequence[-j+1] : gsequence[j-1];

+	g_alt = revp ? gsequencealt[-j+1] : gsequencealt[j-1];

+	if (g == g_alt) {

+	  printf("  %c ",' ');

+	} else {

+	  printf("  %c ",g_alt);

+	}

       }

     }

+    printf("\n");

   }

-  printf("\n");

-#endif

   for (i = 0; i <= rlength; ++i) {

+    printf("%2d ",i);

     if (i == 0) {

       printf("  ");

     } else {

       printf("%c ",revp ? rsequence[-i+1] : rsequence[i-1]);

     }

-    for (j = 0; j <= glength; ++j) {

-      if (j < i - lband) {

-	printf("     ");

-      } else if (j > i + uband) {

-	printf("     ");

-      } else {

-	if (directions_Egap[j][i] == DIAG) {

-	  printf("D");

-	} else {

-	  /* Must be HORIZ */

-	  printf("H");

-	}

-	printf("|");

-	if (directions_nogap[j][i] == DIAG) {

-	  printf("D");

-	} else if (directions_nogap[j][i] == HORIZ) {

-	  printf("H");

+    if (upperp == true) {

+      for (j = 0; j <= glength; ++j) {

+	if (j < i) {

+	  printf("  . ");

+	} else if (j > i + band) {

+	  printf("  . ");

+	} else if (matrix[j][i] < NEG_INFINITY_DISPLAY) {

+	  printf("%3d ",NEG_INFINITY_DISPLAY);

 	} else {

-	  /* Must be VERT */

-	  printf("V");

+	  printf("%3d ",matrix[j][i]);

 	}

-	printf("|");

-	if (directions_Fgap[j][i] == DIAG) {

-	  printf("D");

+      }

+    } else {

+      for (j = 0; j <= glength; ++j) {

+	if (i < j) {

+	  printf("  . ");

+	} else if (i > j + band) {

+	  printf("  . ");

+	} else if (matrix[i][j] < NEG_INFINITY_DISPLAY) {

+	  printf("%3d ",NEG_INFINITY_DISPLAY);

 	} else {

-	  /* Must be VERT */

-	  printf("V");

+	  printf("%3d ",matrix[i][j]);

 	}

       }

-      printf(" ");

     }

     printf("\n");

   }

   printf("\n");

+

   return;

 }

+#endif

+#if defined(DEBUG14) || defined(DEBUG2)

 static void

-Directions16_print (Direction16_T **directions_nogap, Direction16_T **directions_Egap, Direction16_T **directions_Fgap,

-		    int rlength, int glength, char *rsequence, char *gsequence, char *gsequence_alt,

-#ifdef DEBUG14

-		    int goffset, Univcoord_T chroffset, Univcoord_T chrhigh, bool watsonp,

-#endif

-		    bool revp, int lband, int uband) {

+Directions8_print (Direction8_T **directions_nogap, Direction8_T **directions_Egap, Direction8_T **directions_Fgap,

+		   int rlength, int glength, char *rsequence, char *gsequence, char *gsequencealt,

+		   bool revp, int lband, int uband) {

   int i, j;

-  char g_alt;

+  char g, g_alt;

   _mm_lfence();

+  /* j */

+  printf("   ");		/* For i */

+  printf("  ");

+  for (j = 0; j <= glength; ++j) {

+    printf(" %2d   ",j);

+  }

+  printf("\n");

+

   if (gsequence) {

+    printf("   ");		/* For i */

     printf("  ");

     for (j = 0; j <= glength; ++j) {