Hublabel

.gitlab-ci.yml

@@ -158,7 +158,7 @@ oldest-supported-compiler-job:
     GIT_SUBMODULE_STRATEGY: none
     # DO NOT change this version number without updating the README to reflect
     # the requirement bump.
-    COMPILER_VERSION: 9
+    COMPILER_VERSION: 10
 
 
 # We define one job to do the Docker container build

BOTS.md

@@ -0,0 +1,75 @@
+# VG Project Notes
+
+## Building
+- New `.cpp` files auto-discovered
+- Build with `make -j8` or `make obj/whatever.o` to build just one .o.
+- You may be getting errors from `clangd`. If these errors seem spurious, stop and demand a `clangd` that works properly.
+
+## Testing
+
+### Running Bash-TAP Tests
+Use `prove -v` (not `bash`) to execute Bash-TAP tests. This provides proper test harness output and better error reporting.
+
+**Important**: Run `prove` from the `test/` directory:
+```bash
+cd test
+prove -v t/26_deconstruct.t
+```
+
+### Running Unit Tests
+To run all unit tests:
+```bash
+./bin/vg test
+```
+- `./bin/vg test "[tag]"` runs tests matching a tag
+
+#### Writing Unit Tests
+- Framework: Catch v2 (header-only)
+- Include: `#include "catch.hpp"` (in `src/unittest/catch.hpp`)
+- Macros: `TEST_CASE("name", "[tags]")`, `SECTION("name")`, `REQUIRE(cond)`
+- Namespace: `vg::unittest`
+- Directory: `src/unittest/`
+
+### Running All Tests
+```bash
+make test
+```
+
+## Writing Code
+
+### HandleGraph API
+The interfaces in libhandlegraph model a bidirected sequence graph (where nodes have DNA sequences and edges can connect to either the start or end of each involved node).
+
+#### Core types
+- `handle_t` - opaque 64-bit value
+- `nid_t` - node ID type
+- `edge_t` = `pair<handle_t, handle_t>`
+
+#### Key HandleGraph methods
+- `get_handle(nid_t, bool is_reverse=false)` → `handle_t`
+- `get_id(handle_t)` → `nid_t`
+- `get_is_reverse(handle_t)` → `bool`
+- `flip(handle_t)` → `handle_t` (toggle orientation)
+- `get_sequence(handle_t)` → `string` (in handle's orientation)
+- `follow_edges(handle_t, bool go_left, iteratee)` - iterate neighbors
+- `for_each_handle(iteratee, bool parallel=false)` - iterate all nodes
+- `for_each_edge(iteratee, bool parallel=false)` - iterate all edges
+- `has_edge(handle_t left, handle_t right)` → `bool`
+
+#### MutableHandleGraph additions
+- `create_handle(string seq)` / `create_handle(string seq, nid_t id)` → `handle_t`
+- `create_edge(handle_t left, handle_t right)`
+- `destroy_handle(handle_t)` / `destroy_edge(handle_t, handle_t)`
+
+#### HandleGraph algorithms
+- Things like `topological_sort.hpp` and copy_graph.hpp` are in `deps/libhandlegraph/src/include/handlegraph/algorithms`.
+
+#### bdsg::HashGraph
+- Header: `deps/libbdsg/bdsg/include/bdsg/hash_graph.hpp`
+- Implements MutablePathMutableHandleGraph
+- Go-to handlegraph implementation to use
+- In libbdsg
+
+### Utilities
+- `reverse_complement(string)` → `string` in src/utility.hpp
+

CLAUDE.md

@@ -0,0 +1 @@
+BOTS.md

Makefile

@@ -104,7 +104,8 @@ ifeq ($(shell uname -s),Darwin)
     LD_UTIL_RPATH_FLAGS=""
 
     # Homebrew installs a Protobuf that uses an Abseil that is built with C++17, so we need to build with at least C++17
-    CXX_STANDARD?=17
+	# C++20 for spaceship operator and ranges
+    CXX_STANDARD?=20
 
     # We may need libraries from Macports
     ifeq ($(shell if [ -d /opt/local/lib ];then echo 1;else echo 0;fi), 1)
@@ -229,8 +230,9 @@ else
     $(info Compiler $(CXX) is assumed to be GCC)
 
 	# gbwtgraph uses inline variables and our oldest supported compiler has
-	# C++17, so we should use C++17
-    CXX_STANDARD?=17
+	# C++17, so we should use at least C++17.
+    # C++20 for spaceship operator and ranges
+    CXX_STANDARD?=20
 
     # Set an rpath for vg and dependency utils to find installed libraries
     LD_UTIL_RPATH_FLAGS="-Wl,-rpath,$(CWD)/$(LIB_DIR)"
@@ -820,7 +822,7 @@ $(INC_DIR)/dynamic/dynamic.hpp: $(DYNAMIC_DIR)/include/dynamic/*.hpp $(DYNAMIC_D
 	+mkdir -p $(INC_DIR)/dynamic && cp -r $(CWD)/$(DYNAMIC_DIR)/include/dynamic/* $(INC_DIR)/dynamic/
 
 $(INC_DIR)/sparsehash/sparse_hash_map: $(wildcard $(SPARSEHASH_DIR)/**/*.cc) $(wildcard $(SPARSEHASH_DIR)/**/*.h)
-	+cd $(SPARSEHASH_DIR) && ./autogen.sh && LDFLAGS="$(LD_LIB_DIR_FLAGS) $(LDFLAGS)" ./configure --prefix=$(CWD) $(FILTER) && $(MAKE) $(FILTER) && $(MAKE) install
+	+cd $(SPARSEHASH_DIR) && ./autogen.sh && LDFLAGS="$(LD_LIB_DIR_FLAGS) $(LDFLAGS)" ./configure --prefix=$(CWD) $(FILTER) && $(MAKE) src/sparsehash/internal/sparseconfig.h $(FILTER) && $(MAKE) install-data $(FILTER)
 
 $(INC_DIR)/sparsepp/spp.h: $(wildcard $(SPARSEPP_DIR)/sparsepp/*.h)
 	+cp -r $(SPARSEPP_DIR)/sparsepp $(INC_DIR)/

README.md

@@ -93,7 +93,7 @@ On other distros, or if you do not have root access, you will need to perform th
                          liblzma-dev liblz4-dev libffi-dev libcairo-dev libboost-all-dev \
                          libzstd-dev pybind11-dev python3-pybind11 libssl-dev kmc
 
-At present, you will need GCC version 9 or greater, with support for C++17, to compile vg. (Check your version with `gcc --version`.) GCC up to 11.4.0 is supported.
+At present, you will need GCC version 10 or greater, with support for C++20, to compile vg. (Check your version with `gcc --version`.) GCC up to 11.4.0 is supported.
 
 Other libraries may be required. Please report any build difficulties.
 

src/algorithms/chain_items.cpp

@@ -13,6 +13,7 @@
 
 //#define debug_chaining
 //#define debug_transition
+//#define debug_missing_transition
 //#define debug_dp
 
 namespace vg {
@@ -262,6 +263,18 @@
             generate_zip_tree_transitions(seeds, zip_code_tree, max_graph_lookback_bases,
                                           seed_to_starting, seed_to_ending);
 
+#ifdef debug_missing_transition
+        bool has_missing = \
+        find_missing_zip_tree_transitions(seeds, zip_code_tree, max_graph_lookback_bases,
+                                          seed_to_starting, seed_to_ending, distance_index, 
+                                          all_transitions);
+        if (has_missing) {
+            throw std::runtime_error("Zipcode tree iterator failed to output some transitions");
+        } else {
+            cerr << "No missing transitions" << endl;
+        }
+#endif
+
         std::vector<transition_info> filtered_transitions =
             calculate_transition_read_distances(all_transitions, to_chain, max_read_lookback_bases);
 
@@ -385,6 +398,103 @@
     return all_transitions;
 }
 
+bool find_missing_zip_tree_transitions(
+    const std::vector<SnarlDistanceIndexClusterer::Seed>& seeds,
+    const ZipCodeTree& zip_code_tree,
+    size_t max_graph_lookback_bases,
+    const std::unordered_map<size_t, size_t>& seed_to_starting, 
+    const std::unordered_map<size_t, size_t>& seed_to_ending,
+    const SnarlDistanceIndex& distance_index,
+    const std::vector<transition_info>& all_transitions) {
+
+    // {source anchor : {dest anchor : dist}}
+    std::unordered_map<size_t, std::unordered_map<size_t, size_t>> found;
+    for (const auto& transition : all_transitions) {
+        size_t dist_to_save = transition.graph_distance;
+        if (!found.count(transition.from_anchor)) {
+            found[transition.from_anchor] = std::unordered_map<size_t, size_t>();
+        }
+        if (found[transition.from_anchor].count(transition.to_anchor)) {
+            // If a transition appears multiple times, remember the min
+            dist_to_save = std::min(transition.graph_distance, 
+                                    found[transition.from_anchor][transition.to_anchor]);
+        }
+        found[transition.from_anchor][transition.to_anchor] = transition.graph_distance;
+    }
+
+    bool has_missing = false;
+
+    // Helper function to check for a distance between two seeds
+    auto check_distance = [&] (const ZipCodeTree::oriented_seed_t& from_seed, bool rev_from,
+                               const ZipCodeTree::oriented_seed_t& to_seed, bool rev_to) {
+        // XOR to get appropriate orientations
+        rev_from ^= from_seed.is_reversed;
+        rev_to ^= to_seed.is_reversed;
+        if (rev_from != rev_to) {
+            // Cannot be compared; incompatible orientations
+            return;
+        }
+
+        // Look up appropriate anchors
+        auto from_anchor_itr = rev_from ? seed_to_starting.find(from_seed.seed)
+                                        : seed_to_ending.find(from_seed.seed);
+        if ((rev_from && from_anchor_itr == seed_to_starting.end())
+            || (!rev_from && from_anchor_itr == seed_to_ending.end())) {
+            // No anchor exists
+            return;
+        }
+        auto to_anchor_itr = rev_to ? seed_to_ending.find(to_seed.seed)
+                                    : seed_to_starting.find(to_seed.seed);
+        if ((rev_to && to_anchor_itr == seed_to_ending.end())
+            || (!rev_to && to_anchor_itr == seed_to_starting.end())) {
+            // No anchor exists
+            return;
+        }
+
+        // Construct seed positions
+        pos_t from_pos = seeds.at(from_seed.seed).pos;
+        size_t from_length = distance_index.minimum_length(distance_index.get_node_net_handle(id(from_pos)));
+        from_pos = rev_from ? reverse(from_pos, from_length)
+                            : from_pos;
+        pos_t to_pos = seeds.at(to_seed.seed).pos;
+        size_t to_length = distance_index.minimum_length(distance_index.get_node_net_handle(id(to_pos)));
+        to_pos = rev_to ? reverse(to_pos, to_length)
+                        : to_pos;
+
+        // Look up true minimum distance
+        size_t true_distance = minimum_nontrivial_distance(distance_index, from_pos, to_pos);
+        if (true_distance <= max_graph_lookback_bases) {
+            // We should've found this transition
+            auto from_anchor = from_anchor_itr->second;
+            auto to_anchor = to_anchor_itr->second;
+            if (!found.count(from_anchor) 
+                || !found[from_anchor].count(to_anchor)
+                || found[from_anchor][to_anchor] != true_distance) {
+                has_missing = true;
+                cerr << "Missing transition " << from_pos << "->" 
+                     << to_pos << " dist " << true_distance << endl;
+            }
+        }
+    };
+
+    vector<ZipCodeTree::oriented_seed_t> tree_seeds = zip_code_tree.get_all_seeds();
+    for (size_t i = 0; i < tree_seeds.size(); i++) {
+        // Check self-loops
+        check_distance(tree_seeds[i], false, tree_seeds[i], false);
+        check_distance(tree_seeds[i], false, tree_seeds[i], true);
+        check_distance(tree_seeds[i], true, tree_seeds[i], false);
+        for (size_t j = i + 1; j < tree_seeds.size(); j++) {
+            // Check all orientation pairs
+            check_distance(tree_seeds[i], false, tree_seeds[j], false);
+            check_distance(tree_seeds[i], false, tree_seeds[j], true);
+            check_distance(tree_seeds[i], true, tree_seeds[j], false);
+            check_distance(tree_seeds[i], true, tree_seeds[j], true);
+        }
+    }
+
+    return has_missing;
+}
+
 std::vector<transition_info> calculate_transition_read_distances(
     const std::vector<transition_info>& all_transitions,
     const VectorView<Anchor>& to_chain,
@@ -547,6 +657,8 @@
         cerr << "Chaining group of " << to_chain.size() << " items" << endl;
     }
 
+    crash_unless(recomb_penalty >= 0);
+
     // Compute a base seed average length.
     // TODO: Weight anchors differently?
     // TODO: Will this always be the same for all anchors in practice?
@@ -557,6 +669,20 @@
     base_seed_length /= to_chain.size();
 
     chain_scores.resize(to_chain.size());
+
+    // We want to prefer to come from seeds where the transition preserves
+    // access to matching haplotypes, because we don't want to back ourselves
+    // into a corner where we need a recombination when we don't really have
+    // to. So we cheat on the dynamic programming by adding an "evaluation
+    // bonus" to the scores of the different DP options when comparing them. We
+    // keep this bonus out of the actual recorded scores because we don't want
+    // it raising the scores we actually get the more transitions we take.
+    //
+    // We store the bonus used to select the current winning predecessor for
+    // each seed in this vector, which runs alongside the DP table.
+    //
+    // Starting from nowhere means full path conservation, so bonus = recomb_penalty.
+    std::vector<int> eval_bonuses(to_chain.size(), recomb_penalty);
     for (size_t i = 0; i < to_chain.size(); i++) {
         // Set up DP table so we can start anywhere with that item's score, scaled and with bonus applied.
         chain_scores[i] = {(int)(to_chain[i].score() * item_scale + item_bonus), TracedScore::nowhere(), to_chain[i].anchor_end_paths()};
@@ -586,8 +712,20 @@
         }
 
         // If we come from nowhere, we get those points.
-        chain_scores[transition.to_anchor] = std::max(chain_scores[transition.to_anchor], 
-                                                      {(int)item_points, TracedScore::nowhere(), here.anchor_end_paths()});
+        // This also has full path conservation (bonus = recomb_penalty).
+        {
+            TracedScore from_nowhere = {(int)item_points, TracedScore::nowhere(), here.anchor_end_paths()};
+            int nowhere_bonus = recomb_penalty;
+            int eval_nowhere = from_nowhere.score + nowhere_bonus;
+            int eval_current = chain_scores[transition.to_anchor].score + eval_bonuses[transition.to_anchor];
+            if (eval_nowhere > eval_current) {
+                chain_scores[transition.to_anchor] = from_nowhere;
+                eval_bonuses[transition.to_anchor] = nowhere_bonus;
+            } else if (eval_nowhere == eval_current && from_nowhere > chain_scores[transition.to_anchor]) {
+                chain_scores[transition.to_anchor] = from_nowhere;
+                eval_bonuses[transition.to_anchor] = nowhere_bonus;
+            }
+        }
 
         // For each source we could come from
         auto& source = to_chain[transition.from_anchor];
@@ -664,8 +802,34 @@
             TracedScore from_source_score = source_score.add_points(jump_points + item_points)
                                                         .set_shared_paths(here.anchor_paths());
 
-            // Remember that we could make this jump
-            chain_scores[transition.to_anchor] = std::max(chain_scores[transition.to_anchor], from_source_score);
+            // Evaluate heuristic to preserve path flexibility without inflating actual scoring DP.
+            // Bonus = fraction of conserved paths * recomb_penalty.
+            // Bonus is 0 when recombination occurs (no shared paths).
+            int eval_bonus_from = 0;
+            if (recomb_penalty > 0) {
+                int pre_count = __builtin_popcountll(source_score.paths);
+                if (pre_count > 0 && (source_score.paths & here.anchor_start_paths()) != 0) {
+                    // No recombination: bonus = fraction of paths conserved * penalty
+                    int post_count = __builtin_popcountll(from_source_score.paths);
+                    eval_bonus_from = (recomb_penalty * post_count) / pre_count;
+                }
+                // Recombination case (no shared paths): bonus stays 0
+            }
+
+            // Grab the DP table slot we are updating
+            auto& current_best = chain_scores[transition.to_anchor];
+            // Compute the evaluation value for the new candidate
+            int eval_from = from_source_score.score + eval_bonus_from;
+            // Reconstruct the evaluation value for the current winner
+            int eval_best = current_best.score + eval_bonuses[transition.to_anchor];
+
+            if (eval_from > eval_best || (eval_from == eval_best && from_source_score > current_best)) {
+                // Using the evaluation values, and then if tied the real DP
+                // scores, this new candidate beats the previous winner, so
+                // replace it.
+                current_best = from_source_score;
+                eval_bonuses[transition.to_anchor] = eval_bonus_from;
+            }
 
             if (show_work) {
 #ifdef debug_dp

src/algorithms/chain_items.hpp

@@ -477,6 +477,23 @@ std::vector<transition_info> generate_zip_tree_transitions(
     const std::unordered_map<size_t, size_t>& seed_to_starting, 
     const std::unordered_map<size_t, size_t>& seed_to_ending);
 
+/**
+ * Check if all possible transitions were actually found.
+ * 
+ * Iterates over all pairs of seeds and uses the distance index
+ * to determine if there SHOULD have been a transition.
+ * 
+ * Returns if any transitions were missing.
+ */
+bool find_missing_zip_tree_transitions(
+    const std::vector<SnarlDistanceIndexClusterer::Seed>& seeds,
+    const ZipCodeTree& zip_code_tree,
+    size_t max_graph_lookback_bases,
+    const std::unordered_map<size_t, size_t>& seed_to_starting, 
+    const std::unordered_map<size_t, size_t>& seed_to_ending,
+    const SnarlDistanceIndex& distance_index,
+    const std::vector<transition_info>& all_transitions);
+
 /**
  * Calculate read distances for each of the zip tree's transitions.
  * Also filters out transitions that can't be used,