From a7af8a922f8ef5f2e801646903517aa7c32c0868 Mon Sep 17 00:00:00 2001 From: faithokamoto Date: Tue, 18 Aug 2026 15:08:44 -0700 Subject: [PATCH 1/3] as documented --- src/minimizer_mapper_from_chains.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/minimizer_mapper_from_chains.cpp b/src/minimizer_mapper_from_chains.cpp index 543db8095ca..8fccd251c8f 100644 --- a/src/minimizer_mapper_from_chains.cpp +++ b/src/minimizer_mapper_from_chains.cpp @@ -1651,7 +1651,7 @@ void MinimizerMapper::do_chaining_on_trees(const Alignment& aln, const ZipCodeFo get_regular_aligner()->scorer->gap_extension, for_each_transition, scheme, - this->max_alignments, + this->max_chains_per_tree, indel_limit, show_work ); From 8b71fd4b00d7f3c87134ffecdb4b3e1f1c3f76c0 Mon Sep 17 00:00:00 2001 From: faithokamoto Date: Tue, 18 Aug 2026 15:17:24 -0700 Subject: [PATCH 2/3] rename param --- src/minimizer_mapper.hpp | 12 ++++++------ src/minimizer_mapper_from_chains.cpp | 2 +- src/subcommand/giraffe_main.cpp | 16 ++++++++-------- 3 files changed, 15 insertions(+), 15 deletions(-) diff --git a/src/minimizer_mapper.hpp b/src/minimizer_mapper.hpp index 425dcf19ec4..b9ff36610c7 100644 --- a/src/minimizer_mapper.hpp +++ b/src/minimizer_mapper.hpp @@ -305,17 +305,17 @@ class MinimizerMapper : public AlignerClient { static constexpr double default_chain_score_threshold = 100; double chain_score_threshold = default_chain_score_threshold; - /// Disregard the chain score thresholds when they would give us - /// fewer than this many chains aligned. - static constexpr int default_min_chains = 4; - int min_chains = default_min_chains; + /// Disregard the chain_score_threshold (but NOT the min_chain_score) + /// when they would give us fewer than this many chains aligned. + static constexpr int default_target_alignment_attempts = 4; + int target_alignment_attempts = default_target_alignment_attempts; /// Allow up to this many chains per tree static constexpr size_t default_max_chains_per_tree = 1; size_t max_chains_per_tree = default_max_chains_per_tree; - /// Even if we would have fewer than min_chains results, don't - /// process anything with a score smaller than this, per read base. + /// Even if we would have fewer than target_alignment_attempts results, + /// don't process anything with a score smaller than this, per read base. static constexpr double default_min_chain_score_per_base = 0.01; double min_chain_score_per_base = default_min_chain_score_per_base; diff --git a/src/minimizer_mapper_from_chains.cpp b/src/minimizer_mapper_from_chains.cpp index 8fccd251c8f..891556a2d07 100644 --- a/src/minimizer_mapper_from_chains.cpp +++ b/src/minimizer_mapper_from_chains.cpp @@ -1993,7 +1993,7 @@ void MinimizerMapper::do_alignment_on_chains(const Alignment& aln, const std::ve // Go through the chains in estimated-score order. process_until_threshold_b(chain_score_estimates, - chain_score_threshold, min_chains, max_alignments, rng, + chain_score_threshold, target_alignment_attempts, max_alignments, rng, [&](size_t processed_num, size_t item_count) -> bool { // This chain is good enough. // Called in descending score order. diff --git a/src/subcommand/giraffe_main.cpp b/src/subcommand/giraffe_main.cpp index 3dea6124aa9..905acdd7e6a 100644 --- a/src/subcommand/giraffe_main.cpp +++ b/src/subcommand/giraffe_main.cpp @@ -484,10 +484,10 @@ static std::unique_ptr get_options() { double_is_nonnegative ); chaining_opts.add_range( - "min-chains", - &MinimizerMapper::min_chains, - MinimizerMapper::default_min_chains, - "ignore score threshold to get this many chains aligned", + "target-alignment-attempts", + &MinimizerMapper::target_alignment_attempts, + MinimizerMapper::default_target_alignment_attempts, + "ignore score threshold (but not min score) to get this many chains aligned", int_is_nonnegative ); chaining_opts.add_range( @@ -986,7 +986,7 @@ int main_giraffe(int argc, char** argv) { .add_entry("rec-penalty", 2) .add_entry("rec-consistency-bonus", 12) .add_entry("chain-score-threshold", 234.0) - .add_entry("min-chains", 2) + .add_entry("target-alignment-attempts", 2) .add_entry("min-chain-score-per-base", 0.24) .add_entry("max-chains-per-tree", 3) .add_entry("max-min-chain-score", 46) @@ -1045,7 +1045,7 @@ int main_giraffe(int argc, char** argv) { .add_entry("rec-penalty", 2) .add_entry("rec-consistency-bonus", 13) .add_entry("chain-score-threshold", 160.0) - .add_entry("min-chains", 2) + .add_entry("target-alignment-attempts", 2) .add_entry("max-chains-per-tree", 3) // Lowering this can reduce wrong reads in mapping experiments, but // seems to *increase* miscalls; see @@ -1108,7 +1108,7 @@ int main_giraffe(int argc, char** argv) { .add_entry("min-chain-score-per-base", 0.01) .add_entry("max-min-chain-score", 200.0) .add_entry("item-bonus", 0) - .add_entry("min-chains", 3) + .add_entry("target-alignment-attempts", 3) .add_entry("max-chains-per-tree", 5) .add_entry("max-alignments", 4) // Don't use the WFAExtender to connect anchors because it can take tenths of seconds sometimes. @@ -1132,7 +1132,7 @@ int main_giraffe(int argc, char** argv) { .add_entry("max-chaining-problems", std::numeric_limits::max()) .add_entry("max-graph-lookback-bases-per-base", 0) .add_entry("max-indel-bases-per-base", 0) - .add_entry("min-chains", 4) + .add_entry("target-alignment-attempts", 4) .add_entry("max-chains-per-tree", 5) .add_entry("max-alignments", 5); From e67235a28d1df37068032ccddbccf97695f78841 Mon Sep 17 00:00:00 2001 From: faithokamoto Date: Tue, 18 Aug 2026 15:30:28 -0700 Subject: [PATCH 3/3] simplify arguments/returns --- src/algorithms/chain_items.cpp | 163 ++++++++++++--------------- src/algorithms/chain_items.hpp | 44 +++----- src/minimizer_mapper_from_chains.cpp | 2 - src/unittest/chain_items.cpp | 2 +- 4 files changed, 89 insertions(+), 122 deletions(-) diff --git a/src/algorithms/chain_items.cpp b/src/algorithms/chain_items.cpp index aa70cf6e8c5..dcaa7d7bcf8 100644 --- a/src/algorithms/chain_items.cpp +++ b/src/algorithms/chain_items.cpp @@ -382,16 +382,14 @@ int check_recombination(const TracedScore& from, const Anchor& to) { } } -TracedScore chain_items_dp(vector& chain_scores, - const VectorView& to_chain, - const SnarlDistanceIndex& distance_index, - const HandleGraph& graph, - int gap_open, - int gap_extension, - const transition_iterator& for_each_transition, - const ChainScoringScheme& scheme, - size_t max_indel_bases, - bool show_work) { +void chain_items_dp(vector& chain_scores, + const VectorView& to_chain, + const SnarlDistanceIndex& distance_index, + const HandleGraph& graph, + const transition_iterator& for_each_transition, + const ChainScoringScheme& scheme, + size_t max_indel_bases, + bool show_work) { DiagramExplainer diagram(show_work); TSVExplainer dump(show_work, "chaindump"); @@ -590,66 +588,68 @@ TracedScore chain_items_dp(vector& chain_scores, max_indel_bases, iteratee); - - TracedScore best_score = TracedScore::unset(); + if (show_work) { + TracedScore best_score = TracedScore::unset(); - for (size_t to_anchor = 0; to_anchor < to_chain.size(); ++to_anchor) { - // For each destination anchor, now that it is finished, see if it is the winner. - auto& here = to_chain[to_anchor]; + for (size_t to_anchor = 0; to_anchor < to_chain.size(); ++to_anchor) { + // For each destination anchor, now that it is finished, see if it is the winner. + auto& here = to_chain[to_anchor]; - if (show_work) { - cerr << "\tBest way to reach #" << to_anchor << " " << to_chain[to_anchor] - << " is " << chain_scores[to_anchor] << endl; - } - - if (diagram) { - // Draw the item in the diagram - auto item_points = here.score() + scheme.item_bonus; - std::string here_gvnode = "i" + std::to_string(to_anchor); - std::stringstream label_stream; - label_stream << "#" << to_anchor << " " << here << " = " << item_points - << "/" << chain_scores[to_anchor].score; - diagram.add_node(here_gvnode, { - {"label", label_stream.str()} - }); - auto graph_start = here.graph_start(); - std::string graph_gvnode = "n" + std::to_string(id(graph_start)) + (is_rev(graph_start) ? "r" : "f"); - diagram.ensure_node(graph_gvnode, { - {"label", std::to_string(id(graph_start)) + (is_rev(graph_start) ? "-" : "+")}, - {"shape", "box"} - }); - // Show the item as connected to its source graph node - diagram.add_edge(here_gvnode, graph_gvnode, {{"color", "gray"}}); - // Make the next graph node along the same strand - std::string graph_gvnode2 = ("n" + std::to_string(id(graph_start) - + (is_rev(graph_start) ? -1 : 1)) - + (is_rev(graph_start) ? "r" : "f")); - diagram.ensure_node(graph_gvnode2, { - {"label", std::to_string(id(graph_start) + (is_rev(graph_start) ? -1 : 1)) - + (is_rev(graph_start) ? "-" : "+")}, - {"shape", "box"} - }); - // And show them as connected. - diagram.ensure_edge(graph_gvnode, graph_gvnode2, {{"color", "gray"}}); - } - - // See if this is the best overall - best_score.max_in(chain_scores, to_anchor); - - if (show_work) { + if (show_work) { + cerr << "\tBest way to reach #" << to_anchor << " " << to_chain[to_anchor] + << " is " << chain_scores[to_anchor] << endl; + } + + if (diagram) { + // Draw the item in the diagram + auto item_points = here.score() + scheme.item_bonus; + std::string here_gvnode = "i" + std::to_string(to_anchor); + std::stringstream label_stream; + label_stream << "#" << to_anchor << " " << here << " = " << item_points + << "/" << chain_scores[to_anchor].score; + diagram.add_node(here_gvnode, { + {"label", label_stream.str()} + }); + auto graph_start = here.graph_start(); + std::string graph_gvnode = "n" + std::to_string(id(graph_start)) + (is_rev(graph_start) ? "r" : "f"); + diagram.ensure_node(graph_gvnode, { + {"label", std::to_string(id(graph_start)) + (is_rev(graph_start) ? "-" : "+")}, + {"shape", "box"} + }); + // Show the item as connected to its source graph node + diagram.add_edge(here_gvnode, graph_gvnode, {{"color", "gray"}}); + // Make the next graph node along the same strand + std::string graph_gvnode2 = ("n" + std::to_string(id(graph_start) + + (is_rev(graph_start) ? -1 : 1)) + + (is_rev(graph_start) ? "r" : "f")); + diagram.ensure_node(graph_gvnode2, { + {"label", std::to_string(id(graph_start) + (is_rev(graph_start) ? -1 : 1)) + + (is_rev(graph_start) ? "-" : "+")}, + {"shape", "box"} + }); + // And show them as connected. + diagram.ensure_edge(graph_gvnode, graph_gvnode2, {{"color", "gray"}}); + } + + // See if this is the best overall + best_score.max_in(chain_scores, to_anchor); + + if (show_work) { #ifdef debug_dp - cerr << "\tBest chain end so far: " << best_score << endl; + cerr << "\tBest chain end so far: " << best_score << endl; +#endif + } + +#ifdef debug_chaining + std::cerr << "[REC INFO] Recombination number for chain: " << best_score.rec_num + << "\tscore: " << best_score.score << "\tpaths: " << best_score.paths << std::endl; #endif } - } - - return best_score; } vector, int>> chain_items_traceback(const vector& chain_scores, const VectorView& to_chain, - const TracedScore& best_past_ending_score_ever, const ChainScoringScheme& scheme, size_t max_tracebacks) { @@ -678,8 +678,7 @@ vector, int>> chain_items_traceback(const vector traceback; traceback.push_back(trace_from); - // Track the penalty we are off optimal for this traceback - int penalty = best_past_ending_score_ever - chain_scores[trace_from]; + int score = chain_scores[trace_from].score; size_t here = trace_from; #ifdef debug_chaining std::cerr << "[REC INFO] Starting traceback at item #" << here << " with recs: " << chain_scores[here].rec_num << " score: " << chain_scores[here].score << std::endl; @@ -696,9 +695,9 @@ vector, int>> chain_items_traceback(const vector, int>> chain_items_traceback(const vector, int>& a, const std::pair, int>& b) { - // Return true if a has the smaller penalty and belongs first - return a.second < b.second; + // Return true if a has the larger score and belongs first + return a.second > b.second; }); if (tracebacks.size() > max_tracebacks) { @@ -735,8 +734,6 @@ vector, int>> chain_items_traceback(const vector& to_chain, const SnarlDistanceIndex& distance_index, const HandleGraph& graph, - int gap_open, - int gap_extension, const transition_iterator& for_each_transition, const ChainScoringScheme& scheme, size_t max_chains, @@ -755,22 +752,11 @@ ChainsResult find_best_chains(const VectorView& to_chain, // We actually need to do DP vector chain_scores; - TracedScore best_past_ending_score_ever = chain_items_dp(chain_scores, - to_chain, - distance_index, - graph, - gap_open, - gap_extension, - for_each_transition, - scheme, - max_indel_bases, - show_work); -#ifdef debug_chaining - std::cerr << "[REC INFO] Recombination number for chain: " << best_past_ending_score_ever.rec_num << "\tscore: " << best_past_ending_score_ever.score << "\tpaths: " << best_past_ending_score_ever.paths << std::endl; -#endif + chain_items_dp(chain_scores, to_chain, distance_index, graph, + for_each_transition, scheme, max_indel_bases, show_work); // Then do the tracebacks vector, int>> tracebacks = chain_items_traceback( - chain_scores, to_chain, best_past_ending_score_ever, scheme, max_chains); + chain_scores, to_chain, scheme, max_chains); if (tracebacks.empty()) { // Somehow we got nothing @@ -786,8 +772,7 @@ ChainsResult find_best_chains(const VectorView& to_chain, // Everything is already sorted. result.chains.reserve(tracebacks.size()); for (auto& traceback : tracebacks) { - // Move over the list of items and convert penalty to score - int score = best_past_ending_score_ever.score - traceback.second; + // Move over the list of items std::vector chain_indexes = std::move(traceback.first); // Compute the anchor indices in this chain that introduce an @@ -868,7 +853,7 @@ ChainsResult find_best_chains(const VectorView& to_chain, } ChainWithRec entry; - entry.scored_chain = {score, std::move(chain_indexes)}; + entry.scored_chain = {traceback.second, std::move(chain_indexes)}; entry.rec_positions = std::move(rec_positions); entry.rec_intervals = std::move(rec_intervals); result.chains.emplace_back(std::move(entry)); @@ -879,8 +864,6 @@ ChainsResult find_best_chains(const VectorView& to_chain, pair> find_best_chain(const VectorView& to_chain, const SnarlDistanceIndex& distance_index, const HandleGraph& graph, - int gap_open, - int gap_extension, const transition_iterator& for_each_transition, const ChainScoringScheme& scheme, size_t max_indel_bases) { @@ -889,8 +872,6 @@ pair> find_best_chain(const VectorView& to_chain, to_chain, distance_index, graph, - gap_open, - gap_extension, for_each_transition, scheme, 1, diff --git a/src/algorithms/chain_items.hpp b/src/algorithms/chain_items.hpp index 434aa939d6e..209ab0f6c9d 100644 --- a/src/algorithms/chain_items.hpp +++ b/src/algorithms/chain_items.hpp @@ -513,9 +513,7 @@ void add_transition_if_legal(vector& transitions, const VectorV /** * Fill in the given DP table for the explored chain scores ending with each - * item. Returns the best observed score overall from that table, with - * provenance to its location in the table, if tracked in the type. Assumes - * some items exist. + * item. Assumes some items exist. * * We keep all the options to allow us to do multiple tracebacks and find * multiple good (ideally disjoint) chains. @@ -530,36 +528,30 @@ void add_transition_if_legal(vector& transitions, const VectorV * Limits transitions to those involving indels of the given size or less, to * avoid very bad transitions. */ -TracedScore chain_items_dp(vector& chain_scores, - const VectorView& to_chain, - const SnarlDistanceIndex& distance_index, - const HandleGraph& graph, - // TODO: We should maybe just take an EditAlignmentScorer here. - int gap_open, - int gap_extension, - const transition_iterator& for_each_transition, - const ChainScoringScheme& scheme = ChainScoringScheme(), - size_t max_indel_bases = 100, - bool show_work = false - ); +void chain_items_dp(vector& chain_scores, + const VectorView& to_chain, + const SnarlDistanceIndex& distance_index, + const HandleGraph& graph, + const transition_iterator& for_each_transition, + const ChainScoringScheme& scheme = ChainScoringScheme(), + size_t max_indel_bases = 100, + bool show_work = false); /** * Trace back through in the given DP table from the best chain score. * * Returns tracebacks that visit disjoint sets of items, in score order, along - * with their penalties from the optimal score. The best_past_ending_score_ever - * is *not* always the source of the first traceback, if there is a tie. + * with their scores. * - * Tracebacks are constrained to be nonoverlapping by stopping each traceback - * when the optimum place to come from has already been used. The second-best - * place to come from is *not* considered. It might be possible that two - * returned tracebacks could be pasted together to get a higher score, but it - * won't be possible to recombine two tracebacks to get a higher score; no - * edges followed between items will ever need to be cut. + * Tracebacks are constrained to be nonoverlapping by stopping each traceback + * when the optimum place to come from has already been used. The second-best + * place to come from is *not* considered. It might be possible that two + * returned tracebacks could be pasted together to get a higher score, but it + * won't be possible to recombine two tracebacks to get a higher score; no + * edges followed between items will ever need to be cut. */ vector, int>> chain_items_traceback(const vector& chain_scores, const VectorView& to_chain, - const TracedScore& best_past_ending_score_ever, const ChainScoringScheme& scheme = ChainScoringScheme(), size_t max_tracebacks = 1); @@ -576,8 +568,6 @@ vector, int>> chain_items_traceback(const vector& to_chain, const SnarlDistanceIndex& distance_index, const HandleGraph& graph, - int gap_open, - int gap_extension, const transition_iterator& for_each_transition, const ChainScoringScheme& scheme = ChainScoringScheme(), size_t max_chains = 1, @@ -596,8 +586,6 @@ ChainsResult find_best_chains(const VectorView& to_chain, pair> find_best_chain(const VectorView& to_chain, const SnarlDistanceIndex& distance_index, const HandleGraph& graph, - int gap_open, - int gap_extension, const transition_iterator& for_each_transition, const ChainScoringScheme& scheme = ChainScoringScheme(), size_t max_indel_bases = 100); diff --git a/src/minimizer_mapper_from_chains.cpp b/src/minimizer_mapper_from_chains.cpp index 891556a2d07..4ae4a3dc363 100644 --- a/src/minimizer_mapper_from_chains.cpp +++ b/src/minimizer_mapper_from_chains.cpp @@ -1647,8 +1647,6 @@ void MinimizerMapper::do_chaining_on_trees(const Alignment& aln, const ZipCodeFo anchor_view, *distance_index, gbwt_graph, - get_regular_aligner()->scorer->gap_open, - get_regular_aligner()->scorer->gap_extension, for_each_transition, scheme, this->max_chains_per_tree, diff --git a/src/unittest/chain_items.cpp b/src/unittest/chain_items.cpp index cb61c9f4c27..c9a30fc8b3c 100644 --- a/src/unittest/chain_items.cpp +++ b/src/unittest/chain_items.cpp @@ -82,7 +82,7 @@ static pair> run_ziptree_iterator(const HashGraph& graph, // Make iterator for only the first tree // Seriously this is for test cases, only one tree at once - return algorithms::find_best_chain(anchors, distance_index, graph, 6, 1, + return algorithms::find_best_chain(anchors, distance_index, graph, algorithms::zip_tree_transition_iterator(seeds, zip_forest.trees.front(), std::numeric_limits::max(),