Merge remote-tracking branch 'origin/development' into main

tools/align_benchmark/align_benchmark.c

@@ -613,7 +613,9 @@ void align_benchmark_sequential() {
       align_benchmark_print_progress(seqs_processed);
     }
     // DEBUG
-    // mm_allocator_print(stderr,align_input.wf_aligner->mm_allocator,true);
+    //mm_allocator_print(stderr,align_input.wf_aligner->mm_allocator,false);
+    //mm_allocator_print(stderr,align_input.wf_aligner->aligner_forward->mm_allocator,false);
+    //mm_allocator_print(stderr,align_input.wf_aligner->aligner_reverse->mm_allocator,false);
     // Plot
     if (parameters.plot > 0) align_benchmark_plot_wf(&align_input,seqs_processed);
   }
@@ -688,7 +690,7 @@ void align_benchmark_parallel() {
       align_benchmark_print_progress(seqs_processed);
     }
     // DEBUG
-    // mm_allocator_print(stderr,align_input.wf_aligner->mm_allocator,true);
+    //mm_allocator_print(stderr,align_input->wf_aligner->mm_allocator,false);
   }
   // Print benchmark results
   timer_stop(&parameters.timer_global);

wavefront/wavefront_align.c

@@ -70,11 +70,12 @@ bool wavefront_align_reached_limits(
   if (score >= wf_aligner->system.max_alignment_score) {
     wf_aligner->cigar.score = wf_aligner->system.max_alignment_score;
     wf_aligner->align_status.status = WF_STATUS_MAX_SCORE_REACHED;
+    wf_aligner->align_status.score = score;
     return true; // Stop
   }
   // Global probing interval
   alignment_system_t* const system = &wf_aligner->system;
-  if ((score%system->probe_interval_global) != 0) return false; // Continue
+  if (score % system->probe_interval_global != 0) return false; // Continue
   if (system->verbose) {
     wavefront_aligner_print_status(stderr,wf_aligner,score); // DEBUG
   }
@@ -102,6 +103,7 @@ bool wavefront_align_reached_limits(
   const uint64_t wf_memory_used = wavefront_aligner_get_size(wf_aligner);
   if (wf_memory_used > system->max_memory_abort) {
     wf_aligner->align_status.status = WF_STATUS_OOM;
+    wf_aligner->align_status.score = score;
     return true; // Stop
   }
   // Otherwise continue
@@ -290,6 +292,8 @@ void wavefront_align_terminate(
     wf_aligner->cigar.score = (distance_metric <= edit) ? score :
         WF_PENALTIES_GET_SW_SCORE(swg_match_score,pattern_length,text_length,score);
   }
+  // Set successful
+  wf_aligner->align_status.status = WF_STATUS_SUCCESSFUL;
 }
 /*
  * General Alignment
@@ -309,20 +313,16 @@ int wavefront_align_sequences(
     if (finished) {
       // DEBUG
       // wavefront_aligner_print(stderr,wf_aligner,0,score,7,0);
-      if (align_status->status == WF_STATUS_SUCCESSFUL) {
+      if (align_status->status == WF_STATUS_END_REACHED) {
         wavefront_align_terminate(wf_aligner,score);
       }
-      align_status->score = score;
       return align_status->status;
     }
     // Compute (s+1)-wavefront
     ++score;
     (*wf_align_compute)(wf_aligner,score);
     // Probe limits
-    if (wavefront_align_reached_limits(wf_aligner,score)) {
-      align_status->score = score;
-      return align_status->status;
-    }
+    if (wavefront_align_reached_limits(wf_aligner,score)) return align_status->status;
     // PROFILE
     if (wf_aligner->plot_params.plot_enabled) {
       wavefront_plot(wf_aligner,wf_aligner->pattern,wf_aligner->text,score);

wavefront/wavefront_aligner.c

@@ -369,38 +369,71 @@ void wavefront_aligner_set_alignment_free_ends(
 void wavefront_aligner_set_heuristic_none(
     wavefront_aligner_t* const wf_aligner) {
   wavefront_heuristic_set_none(&wf_aligner->heuristic);
+  if (wf_aligner->bidirectional_alignment) {
+    wavefront_heuristic_set_none(&wf_aligner->aligner_forward->heuristic);
+    wavefront_heuristic_set_none(&wf_aligner->aligner_reverse->heuristic);
+  }
 }
 void wavefront_aligner_set_heuristic_banded_static(
     wavefront_aligner_t* const wf_aligner,
     const int band_min_k,
     const int band_max_k) {
   wavefront_heuristic_set_banded_static(&wf_aligner->heuristic,band_min_k,band_max_k);
+  if (wf_aligner->bidirectional_alignment) {
+    wavefront_heuristic_set_banded_static(&wf_aligner->aligner_forward->heuristic,band_min_k,band_max_k);
+    wavefront_heuristic_set_banded_static(&wf_aligner->aligner_reverse->heuristic,band_min_k,band_max_k);
+  }
 }
 void wavefront_aligner_set_heuristic_banded_adaptive(
     wavefront_aligner_t* const wf_aligner,
     const int band_min_k,
     const int band_max_k,
     const int score_steps) {
-  wavefront_heuristic_set_banded_adaptive(&wf_aligner->heuristic,band_min_k,band_max_k,score_steps);
+  wavefront_heuristic_set_banded_adaptive(
+      &wf_aligner->heuristic,band_min_k,band_max_k,score_steps);
+  if (wf_aligner->bidirectional_alignment) {
+    wavefront_heuristic_set_banded_adaptive(
+        &wf_aligner->aligner_forward->heuristic,band_min_k,band_max_k,score_steps);
+    wavefront_heuristic_set_banded_adaptive(
+        &wf_aligner->aligner_reverse->heuristic,band_min_k,band_max_k,score_steps);
+  }
 }
 void wavefront_aligner_set_heuristic_wfadaptive(
     wavefront_aligner_t* const wf_aligner,
     const int min_wavefront_length,
     const int max_distance_threshold,
     const int score_steps) {
-  wavefront_heuristic_set_wfadaptive(&wf_aligner->heuristic,min_wavefront_length,max_distance_threshold,score_steps);
+  wavefront_heuristic_set_wfadaptive(
+      &wf_aligner->heuristic,
+      min_wavefront_length,max_distance_threshold,score_steps);
+  if (wf_aligner->bidirectional_alignment) {
+    wavefront_heuristic_set_wfadaptive(
+        &wf_aligner->aligner_forward->heuristic,
+        min_wavefront_length,max_distance_threshold,score_steps);
+    wavefront_heuristic_set_wfadaptive(
+        &wf_aligner->aligner_reverse->heuristic,
+        min_wavefront_length,max_distance_threshold,score_steps);
+  }
 }
 void wavefront_aligner_set_heuristic_xdrop(
     wavefront_aligner_t* const wf_aligner,
     const int xdrop,
     const int score_steps) {
   wavefront_heuristic_set_xdrop(&wf_aligner->heuristic,xdrop,score_steps);
+  if (wf_aligner->bidirectional_alignment) {
+    wavefront_heuristic_set_xdrop(&wf_aligner->aligner_forward->heuristic,xdrop,score_steps);
+    wavefront_heuristic_set_xdrop(&wf_aligner->aligner_reverse->heuristic,xdrop,score_steps);
+  }
 }
 void wavefront_aligner_set_heuristic_zdrop(
     wavefront_aligner_t* const wf_aligner,
     const int ydrop,
     const int score_steps) {
   wavefront_heuristic_set_zdrop(&wf_aligner->heuristic,ydrop,score_steps);
+  if (wf_aligner->bidirectional_alignment) {
+    wavefront_heuristic_set_zdrop(&wf_aligner->aligner_forward->heuristic,ydrop,score_steps);
+    wavefront_heuristic_set_zdrop(&wf_aligner->aligner_reverse->heuristic,ydrop,score_steps);
+  }
 }
 /*
  * Match-funct configuration

wavefront/wavefront_aligner.h

@@ -52,6 +52,7 @@
  */
 #define WF_STATUS_SUCCESSFUL               0
 #define WF_STATUS_IN_PROGRESS              1
+#define WF_STATUS_END_REACHED              2 /* Internal */
 #define WF_STATUS_UNFEASIBLE              -1
 #define WF_STATUS_MAX_SCORE_REACHED       -2
 #define WF_STATUS_OOM                     -3

wavefront/wavefront_bialign.c

@@ -381,9 +381,9 @@ int wavefront_bialign_find_breakpoint(
   // Prepare and perform first bialignment step
   breakpoint->score = INT_MAX;
   end_reached = wavefront_extend_end2end_max(wf_aligner_forward,score_forward,&forward_max_ak);
-  if (end_reached) return WF_STATUS_UNFEASIBLE;
+  if (end_reached) return wf_aligner_forward->align_status.status;
   end_reached = wavefront_extend_end2end_max(wf_aligner_reverse,score_reverse,&reverse_max_ak);
-  if (end_reached) return WF_STATUS_UNFEASIBLE;
+  if (end_reached) return wf_aligner_reverse->align_status.status;
   // Compute wavefronts of increasing score until both wavefronts overlap
   int max_ak = 0;
   bool last_wf_forward;
@@ -394,7 +394,7 @@ int wavefront_bialign_find_breakpoint(
     ++score_forward;
     (*wf_align_compute)(wf_aligner_forward,score_forward);
     end_reached = wavefront_extend_end2end_max(wf_aligner_forward,score_forward,&max_ak);
-    if (end_reached) return WF_STATUS_UNFEASIBLE;
+    if (end_reached) return wf_aligner_forward->align_status.status;
     if (forward_max_ak < max_ak) forward_max_ak = max_ak;
     last_wf_forward = true;
     // Check if they are close to collision
@@ -403,7 +403,7 @@ int wavefront_bialign_find_breakpoint(
     ++score_reverse;
     (*wf_align_compute)(wf_aligner_reverse,score_reverse);
     end_reached = wavefront_extend_end2end_max(wf_aligner_reverse,score_reverse,&max_ak);
-    if (end_reached) return WF_STATUS_UNFEASIBLE;
+    if (end_reached) return wf_aligner_reverse->align_status.status;
     if (reverse_max_ak < max_ak) reverse_max_ak = max_ak;
     last_wf_forward = false;
     // Check max-score reached (to stop)
@@ -422,7 +422,7 @@ int wavefront_bialign_find_breakpoint(
       ++score_reverse;
       (*wf_align_compute)(wf_aligner_reverse,score_reverse);
       end_reached = wavefront_extend_end2end(wf_aligner_reverse,score_reverse);
-      if (end_reached) return WF_STATUS_UNFEASIBLE;
+      if (end_reached) return wf_aligner_reverse->align_status.status;
     }
     // Check overlapping wavefronts
     const int min_score_forward = (score_forward > max_score_scope-1) ? score_forward - (max_score_scope-1) : 0;
@@ -432,7 +432,7 @@ int wavefront_bialign_find_breakpoint(
     ++score_forward;
     (*wf_align_compute)(wf_aligner_forward,score_forward);
     end_reached = wavefront_extend_end2end(wf_aligner_forward,score_forward);
-    if (end_reached) return WF_STATUS_UNFEASIBLE;
+    if (end_reached) return wf_aligner_forward->align_status.status;
     // Check max-score reached (to stop)
     if (score_reverse + score_forward >= max_alignment_score) return WF_STATUS_MAX_SCORE_REACHED;
     // Enable always
@@ -499,6 +499,48 @@ void wavefront_bialign_base(
       pattern,pattern_length,
       text,text_length,true);
   cigar_append(cigar,&wf_aligner->cigar);
+  if (rlevel == 0) cigar->score = wf_aligner->cigar.score;
+}
+void wavefront_bialign_exception(
+    wavefront_aligner_t* const wf_aligner,
+    const char* const pattern,
+    const int pattern_length,
+    const char* const text,
+    const int text_length,
+    alignment_form_t* const form,
+    const affine2p_matrix_type component_begin,
+    const affine2p_matrix_type component_end,
+    cigar_t* const cigar,
+    const int rlevel,
+    const int align_status) {
+  // Check max-score reached or unfeasible alignment
+  if (align_status == WF_STATUS_MAX_SCORE_REACHED ||
+      align_status == WF_STATUS_UNFEASIBLE) {
+    wf_aligner->align_status.status = align_status;
+    return;
+  }
+  // Check end reached
+  if (align_status == WF_STATUS_END_REACHED) {
+    // Retrieve score when end was reached
+    int score_reached;
+    if (wf_aligner->aligner_forward->align_status.status == WF_STATUS_END_REACHED) {
+      score_reached = wf_aligner->aligner_forward->align_status.score;
+    } else {
+      score_reached = wf_aligner->aligner_reverse->align_status.score;
+    }
+    // Fallback if possible
+    if (score_reached <= WF_BIALIGN_FALLBACK_MIN_SCORE) {
+      wavefront_bialign_base(
+          wf_aligner,pattern,pattern_length,text,text_length,
+          form,component_begin,component_end,cigar,rlevel);
+    } else {
+      wf_aligner->align_status.status = WF_STATUS_UNFEASIBLE;
+    }
+    return;
+  }
+  // Otherwise
+  fprintf(stderr,"[WFA::BiAlign] Unknown condition \n");
+  exit(1);
 }
 void wavefront_bialign_alignment(
     wavefront_aligner_t* const wf_aligner,
@@ -534,16 +576,13 @@ void wavefront_bialign_alignment(
       pattern,pattern_length,text,text_length,
       wf_aligner->penalties.distance_metric,
       form,component_begin,component_end,&breakpoint);
-  if (align_status == WF_STATUS_UNFEASIBLE) {
-    wavefront_bialign_base(
+  if (align_status != WF_STATUS_SUCCESSFUL) {
+    wavefront_bialign_exception(
         wf_aligner,pattern,pattern_length,text,text_length,
-        form,component_begin,component_end,cigar,rlevel);
-    return;
-  }
-  if (align_status == WF_STATUS_MAX_SCORE_REACHED) {
-    wf_aligner->align_status.status = WF_STATUS_MAX_SCORE_REACHED;
+        form,component_begin,component_end,cigar,rlevel,align_status);
     return;
   }
+  // Breakpoint found
   const int breakpoint_h = WAVEFRONT_H(breakpoint.k_forward,breakpoint.offset_forward);
   const int breakpoint_v = WAVEFRONT_V(breakpoint.k_forward,breakpoint.offset_forward);
   // DEBUG
@@ -587,20 +626,24 @@ void wavefront_bialign_compute_score(
       pattern,pattern_length,text,text_length,
       wf_aligner->penalties.distance_metric,
       &wf_aligner->alignment_form,affine_matrix_M,affine_matrix_M,&breakpoint);
-  if (align_status == WF_STATUS_UNFEASIBLE) {
-    wf_aligner->align_status.status = WF_STATUS_UNFEASIBLE;
+  if (align_status == WF_STATUS_MAX_SCORE_REACHED || align_status == WF_STATUS_UNFEASIBLE) {
+    wf_aligner->align_status.status = align_status;
     return;
   }
-  if (align_status == WF_STATUS_MAX_SCORE_REACHED) {
-    wf_aligner->align_status.status = WF_STATUS_MAX_SCORE_REACHED;
-    return;
+  if (align_status == WF_STATUS_END_REACHED) {
+    if (wf_aligner->aligner_forward->align_status.status == WF_STATUS_END_REACHED) {
+      breakpoint.score = wf_aligner->aligner_forward->align_status.score;
+    } else {
+      breakpoint.score = wf_aligner->aligner_reverse->align_status.score;
+    }
   }
   // Report score
   const distance_metric_t distance_metric = wf_aligner->penalties.distance_metric;
   const int swg_match_score = wf_aligner->penalties.match;
   cigar_clear(&wf_aligner->cigar);
   wf_aligner->cigar.score = (distance_metric <= edit) ? breakpoint.score :
       WF_PENALTIES_GET_SW_SCORE(swg_match_score,pattern_length,text_length,breakpoint.score);
+  wf_aligner->align_status.status = WF_STATUS_SUCCESSFUL;
 }
 /*
  * Bidirectional dispatcher

wavefront/wavefront_extend.c

@@ -309,6 +309,7 @@ int wavefront_extend_end2end_max(
     // Check alignment feasibility (for heuristic variants that can lead to no solution)
     if (wf_aligner->align_status.num_null_steps > wf_aligner->wf_components.max_score_scope) {
       wf_aligner->align_status.status = WF_STATUS_UNFEASIBLE;
+      wf_aligner->align_status.score = score;
       return 1; // Done
     }
     return 0; // Not done
@@ -342,7 +343,8 @@ int wavefront_extend_end2end_max(
   // Check end-to-end finished
   const bool end_reached = wavefront_extend_end2end_check_termination(wf_aligner,mwavefront,score,score_mod);
   if (end_reached) {
-    wf_aligner->align_status.status = WF_STATUS_SUCCESSFUL;
+    wf_aligner->align_status.status = WF_STATUS_END_REACHED;
+    wf_aligner->align_status.score = score;
     return 1; // Done
   }
   // Cut-off wavefront heuristically
@@ -365,6 +367,7 @@ int wavefront_extend_end2end(
     // Check alignment feasibility (for heuristic variants that can lead to no solution)
     if (wf_aligner->align_status.num_null_steps > wf_aligner->wf_components.max_score_scope) {
       wf_aligner->align_status.status = WF_STATUS_UNFEASIBLE;
+      wf_aligner->align_status.score = score;
       return 1; // Done
     }
     return 0; // Not done
@@ -392,7 +395,8 @@ int wavefront_extend_end2end(
   // Check end-to-end finished
   end_reached = wavefront_extend_end2end_check_termination(wf_aligner,mwavefront,score,score_mod);
   if (end_reached) {
-    wf_aligner->align_status.status = WF_STATUS_SUCCESSFUL;
+    wf_aligner->align_status.status = WF_STATUS_END_REACHED;
+    wf_aligner->align_status.score = score;
     return 1; // Done
   }
   // Cut-off wavefront heuristically
@@ -414,6 +418,7 @@ int wavefront_extend_endsfree(
     // Check alignment feasibility (for heuristic variants that can lead to no solution)
     if (wf_aligner->align_status.num_null_steps > wf_aligner->wf_components.max_score_scope) {
       wf_aligner->align_status.status = WF_STATUS_UNFEASIBLE;
+      wf_aligner->align_status.score = score;
       return 1; // Done
     }
     return 0; // Not done
@@ -441,7 +446,8 @@ int wavefront_extend_endsfree(
 #endif
   }
   if (end_reached) {
-    wf_aligner->align_status.status = WF_STATUS_SUCCESSFUL;
+    wf_aligner->align_status.status = WF_STATUS_END_REACHED;
+    wf_aligner->align_status.score = score;
     return 1; // Done
   }
   // Cut-off wavefront heuristically
@@ -463,6 +469,7 @@ int wavefront_extend_custom(
     // Check alignment feasibility (for heuristic variants that can lead to no solution)
     if (wf_aligner->align_status.num_null_steps > wf_aligner->wf_components.max_score_scope) {
       wf_aligner->align_status.status = WF_STATUS_UNFEASIBLE;
+      wf_aligner->align_status.score = score;
       return 1; // Done
     }
     return 0; // Not done
@@ -495,7 +502,8 @@ int wavefront_extend_custom(
     end_reached = wavefront_extend_end2end_check_termination(wf_aligner,mwavefront,score,score_mod);
   }
   if (end_reached) {
-    wf_aligner->align_status.status = WF_STATUS_SUCCESSFUL;
+    wf_aligner->align_status.status = WF_STATUS_END_REACHED;
+    wf_aligner->align_status.score = score;
     return 1; // Done
   }
   // Cut-off wavefront heuristically