Making room for BiWFA (support functions)

Makefile

@@ -44,7 +44,7 @@ ifeq ($(BUILD_EXAMPLES),1)
     APPS+=examples
 endif
 
-all: CC_FLAGS+=-O3 -march=native -flto
+all: CC_FLAGS+=-O3 -march=native #-flto -ffat-lto-objects
 all: build
 
 debug: build

alignment/cigar.c

@@ -251,6 +251,41 @@ void cigar_copy(
          cigar_src->operations+cigar_src->begin_offset,
          cigar_src->end_offset-cigar_src->begin_offset);
 }
+void cigar_append(
+    cigar_t* const cigar_dst,
+    cigar_t* const cigar_src) {
+  // Append
+  const int cigar_length = cigar_src->end_offset - cigar_src->begin_offset;
+  char* const operations_src = cigar_src->operations + cigar_src->begin_offset;
+  char* const operations_dst = cigar_dst->operations + cigar_dst->end_offset;
+  memcpy(operations_dst,operations_src,cigar_length);
+  // Update offset
+  cigar_dst->end_offset += cigar_length;
+}
+void cigar_append_deletion(
+    cigar_t* const cigar,
+    const int length) {
+  // Append deletions
+  char* const operations = cigar->operations + cigar->end_offset;
+  int i;
+  for (i=0;i<length;++i) {
+    operations[i] = 'D';
+  }
+  // Update offset
+  cigar->end_offset += length;
+}
+void cigar_append_insertion(
+    cigar_t* const cigar,
+    const int length) {
+  // Append insertions
+  char* const operations = cigar->operations + cigar->end_offset;
+  int i;
+  for (i=0;i<length;++i) {
+    operations[i] = 'I';
+  }
+  // Update offset
+  cigar->end_offset += length;
+}
 bool cigar_check_alignment(
     FILE* const stream,
     const char* const pattern,

alignment/cigar.h

@@ -104,6 +104,17 @@ int cigar_cmp(
 void cigar_copy(
     cigar_t* const cigar_dst,
     cigar_t* const cigar_src);
+
+void cigar_append(
+    cigar_t* const cigar_dst,
+    cigar_t* const cigar_src);
+void cigar_append_deletion(
+    cigar_t* const cigar,
+    const int length);
+void cigar_append_insertion(
+    cigar_t* const cigar,
+    const int length);
+
 bool cigar_check_alignment(
     FILE* const stream,
     const char* const pattern,

scripts/wfa.alg.cmp.py

@@ -117,8 +117,9 @@ def plot_score_distribution(stats,input_path1,input_path2):
   #   ax1.xaxis.grid(True)
   #   ax1.yaxis.grid(True)
   # Plot score histogram
-  ax1.hist(stats.scores1,color="royalblue",edgecolor = 'black',alpha=0.75)
-  ax1.hist(stats.scores2,color="darkorange",edgecolor = 'black',alpha=0.75)
+  b=range(min(stats.scores1+stats.scores2), max(stats.scores1+stats.scores2) + 10, 10)
+  ax1.hist(stats.scores1,bins=b,color="royalblue",edgecolor = 'black',alpha=0.75)
+  ax1.hist(stats.scores2,bins=b,color="darkorange",edgecolor = 'black',alpha=0.75)
   # Leyend
   handles = [Rectangle((0,0),1,1,color=c,ec="k") for c in ["royalblue","darkorange"]]
   labels= [input_path1,input_path2]
@@ -155,11 +156,11 @@ def compare_alignments(input_path1,input_path2,penalties,use_score,ignore_misms,
       if use_score:
         cigar1 = None
         cigar2 = None
-        score1 = fields1[0] if len(fields1)==2 else fields1[2]
-        score2 = fields2[0] if len(fields2)==2 else fields2[2]
+        score1 = fields1[0] if len(fields1)<=2 else fields1[2]
+        score2 = fields2[0] if len(fields2)<=2 else fields2[2]
       else:
-        cigar1 = fields1[1] if len(fields1)==2 else fields1[5]
-        cigar2 = fields2[1] if len(fields2)==2 else fields2[5]
+        cigar1 = fields1[1] if len(fields1)<=2 else fields1[5]
+        cigar2 = fields2[1] if len(fields2)<=2 else fields2[5]
         # Evaluate CIGAR's score
         score1 = cigar_compute_score(cigar1,penalties,ignore_misms)
         score2 = cigar_compute_score(cigar2,penalties,ignore_misms)

scripts/wfa.utest.cmp.alignments.sh

@@ -17,15 +17,21 @@ do
   FILE_ALG2="$FOLDER2/$FILENAME"
   echo -ne "[UTest::$PREFIX]  \t"
   # Check existence
-  if [[ ! -f "$FILE_ALG2" ]]; then
+  if [[ ! -f "$FILE_ALG2" ]]
+  then
     echo "$FILE_ALG2 doesn't exist."
     continue
   fi
   # Check diff
   if [[ $(diff $FILE_ALG1 $FILE_ALG2) ]] 
   then
-    echo "Error"
-    continue
+    if [[ $(diff  <(awk '{print $1}' $FILE_ALG1) <(awk '{print $1}' $FILE_ALG2)) ]]
+    then
+      echo "Error"
+      continue
+    else
+      echo -n "OK-Score"
+    fi
   else
     echo -n "OK"
   fi

tools/align_benchmark/align_benchmark.c

@@ -798,7 +798,7 @@ void parse_arguments(int argc,char** argv) {
     /* Other alignment parameters */
     { "bandwidth", required_argument, 0, 2000 },
     /* Misc */
-    { "check", optional_argument, 0, 'c' },
+    { "check", required_argument, 0, 'c' },
     { "check-distance", required_argument, 0, 3001 },
     { "check-bandwidth", required_argument, 0, 3002 },
     { "plot", optional_argument, 0, 3003 },
@@ -996,9 +996,9 @@ void parse_arguments(int argc,char** argv) {
      * Misc
      */
     case 'c':
-      if (optarg ==  NULL) { // default = score
+      if (optarg ==  NULL) { // default = correct
         parameters.check_correct = true;
-        parameters.check_score = true;
+        parameters.check_score = false;
         parameters.check_alignments = false;
       } else if (strcasecmp(optarg,"display")==0) {
         parameters.check_display = true;

tools/align_benchmark/benchmark/benchmark_utils.c

@@ -107,7 +107,7 @@ void benchmark_print_output_lite(
     const int score,
     cigar_t* const cigar) {
   // Retrieve CIGAR
-  char* const cigar_str = malloc(cigar->end_offset-cigar->begin_offset);
+  char* const cigar_str = malloc(2*(cigar->end_offset-cigar->begin_offset));
   cigar_sprint(cigar_str,cigar,true);
   // Print
   fprintf(stream,"%d\t%s\n",score,cigar_str);
@@ -120,7 +120,7 @@ void benchmark_print_output_full(
     const int score,
     cigar_t* const cigar) {
   // Retrieve CIGAR
-  char* const cigar_str = malloc(cigar->end_offset-cigar->begin_offset);
+  char* const cigar_str = malloc(2*(cigar->end_offset-cigar->begin_offset));
   cigar_sprint(cigar_str,cigar,true);
   // Print
   fprintf(stream,"%d\t%d\t%d\t%s\t%s\t%s\n",

utils/string_padded.c

@@ -43,6 +43,7 @@ void strings_padded_add_padding(
     const char padding_value,
     char** const buffer_padded,
     char** const buffer_padded_begin,
+    const bool reverse_sequence,
     mm_allocator_t* const mm_allocator) {
   // Allocate
   const int buffer_padded_length = begin_padding_length + buffer_length + end_padding_length;
@@ -51,7 +52,14 @@ void strings_padded_add_padding(
   memset(*buffer_padded,padding_value,begin_padding_length);
   // Copy buffer
   *buffer_padded_begin = *buffer_padded + begin_padding_length;
-  memcpy(*buffer_padded_begin,buffer,buffer_length);
+  if (reverse_sequence) {
+    int i;
+    for (i=0;i<buffer_length;i++) {
+      (*buffer_padded_begin)[i] = buffer[buffer_length-1-i];
+    }
+  } else {
+    memcpy(*buffer_padded_begin,buffer,buffer_length);
+  }
   // Add end padding
   memset(*buffer_padded_begin+buffer_length,padding_value,end_padding_length);
 }
@@ -61,6 +69,7 @@ strings_padded_t* strings_padded_new(
     const char* const text,
     const int text_length,
     const int padding_length,
+    const bool reverse_sequences,
     mm_allocator_t* const mm_allocator) {
   // Allocate
   strings_padded_t* const strings_padded =
@@ -76,12 +85,14 @@ strings_padded_t* strings_padded_new(
       pattern,pattern_length,
       pattern_begin_padding_length,pattern_end_padding_length,'?',
       &(strings_padded->pattern_padded_buffer),
-      &(strings_padded->pattern_padded),mm_allocator);
+      &(strings_padded->pattern_padded),
+      reverse_sequences,mm_allocator);
   strings_padded_add_padding(
       text,text_length,
       text_begin_padding_length,text_end_padding_length,'!',
       &(strings_padded->text_padded_buffer),
-      &(strings_padded->text_padded),mm_allocator);
+      &(strings_padded->text_padded),
+      reverse_sequences,mm_allocator);
   // Return
   return strings_padded;
 }
@@ -91,6 +102,7 @@ strings_padded_t* strings_padded_new_rhomb(
     const char* const text,
     const int text_length,
     const int padding_length,
+    const bool reverse_sequences,
     mm_allocator_t* const mm_allocator) {
   // Allocate
   strings_padded_t* const strings_padded =
@@ -106,12 +118,14 @@ strings_padded_t* strings_padded_new_rhomb(
       pattern,pattern_length,
       pattern_begin_padding_length,pattern_end_padding_length,'?',
       &(strings_padded->pattern_padded_buffer),
-      &(strings_padded->pattern_padded),mm_allocator);
+      &(strings_padded->pattern_padded),
+      reverse_sequences,mm_allocator);
   strings_padded_add_padding(
       text,text_length,
       text_begin_padding_length,text_end_padding_length,'!',
       &(strings_padded->text_padded_buffer),
-      &(strings_padded->text_padded),mm_allocator);
+      &(strings_padded->text_padded),
+      reverse_sequences,mm_allocator);
   // Set lengths
   strings_padded->pattern_length = pattern_length;
   strings_padded->text_length = text_length;

utils/string_padded.h

@@ -63,13 +63,15 @@ strings_padded_t* strings_padded_new(
     const char* const text,
     const int text_length,
     const int padding_length,
+    const bool reverse_sequences,
     mm_allocator_t* const mm_allocator);
 strings_padded_t* strings_padded_new_rhomb(
     const char* const pattern,
     const int pattern_length,
     const char* const text,
     const int text_length,
     const int padding_length,
+    const bool reverse_sequences,
     mm_allocator_t* const mm_allocator);
 void strings_padded_delete(
     strings_padded_t* const strings_padded);

wavefront/wavefront.c

@@ -92,7 +92,6 @@ void wavefront_init(
   wavefront->null = false;
   wavefront->lo =  1;
   wavefront->hi = -1;
-  wavefront->k_alignment_end = WAVEFRONT_DIAGONAL_NULL;
   wavefront->bt_occupancy_max = 0;
   // Setup elements
   wavefront->offsets = wavefront->offsets_mem - min_lo; // Center at k=0
@@ -114,7 +113,6 @@ void wavefront_init_null(
   wavefront->null = true;
   wavefront->lo =  1;
   wavefront->hi = -1;
-  wavefront->k_alignment_end = WAVEFRONT_DIAGONAL_NULL;
   wavefront->bt_occupancy_max = 0;
   // Setup elements
   wavefront->offsets = wavefront->offsets_mem - min_lo; // Center at k=0

wavefront/wavefront.h

@@ -37,6 +37,15 @@
 #include "wavefront_offset.h"
 #include "wavefront_backtrace_buffer.h"
 
+/*
+ * Alignment position
+ */
+typedef struct {
+  int score;          // Score
+  int k;              // Diagonal
+  wf_offset_t offset; // Offset
+} wavefront_pos_t;
+
 /*
  * Wavefront
  */
@@ -50,7 +59,6 @@ typedef struct {
   bool null;                           // Is null interval?
   int lo;                              // Lowest diagonal (inclusive)
   int hi;                              // Highest diagonal (inclusive)
-  int k_alignment_end;                 // Mark WF's diagonal that reached the end of the alignment (semi-global)
   int bt_occupancy_max;                // Maximum number of pcigar-ops stored on the Backtrace-block
   // Wavefront elements
   wf_offset_t* offsets;                // Offsets (k-centered)