diff --git a/tools/align_benchmark/align_benchmark.c b/tools/align_benchmark/align_benchmark.c index 1025ed0..50b36ba 100644 --- a/tools/align_benchmark/align_benchmark.c +++ b/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(¶meters.timer_global); diff --git a/wavefront/wavefront_align.c b/wavefront/wavefront_align.c index 0199d8a..4f88e62 100644 --- a/wavefront/wavefront_align.c +++ b/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); diff --git a/wavefront/wavefront_aligner.c b/wavefront/wavefront_aligner.c index 9ba6160..a3356bb 100644 --- a/wavefront/wavefront_aligner.c +++ b/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 diff --git a/wavefront/wavefront_aligner.h b/wavefront/wavefront_aligner.h index bce8028..8c396f6 100644 --- a/wavefront/wavefront_aligner.h +++ b/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 diff --git a/wavefront/wavefront_bialign.c b/wavefront/wavefront_bialign.c index e67e5f7..f542bcc 100644 --- a/wavefront/wavefront_bialign.c +++ b/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,13 +626,16 @@ 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; @@ -601,6 +643,7 @@ void wavefront_bialign_compute_score( 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 diff --git a/wavefront/wavefront_extend.c b/wavefront/wavefront_extend.c index 6f82f88..1b83924 100644 --- a/wavefront/wavefront_extend.c +++ b/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