From 8935f4e453d6beb5eae4771eacec4581c1419795 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sun, 1 Nov 2026 23:03:19 +0100 Subject: [PATCH 01/52] Add the traversal label oracle and /resolve support profiling MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit First piece of the labeled-traversal v1 stack, rebuilt from Gunnar's gr/labeled-traversal (d5d26baa) without the budget, pacing, trace and annotate layers: - traversal_types: regimes, label references, arms and end reasons. - LabelOracle: the oriented graph to walk (a per-request CanonicalDBG in the PRIMARY regime), node-to-annotation-key mapping through AnnotationKeyMapper, label resolution (column, then header), and batched get_rows / get_row_tuples. - LabelQuery: membership of a fixed permitted label set per key, column labels from rows, header labels from tuple rows by coordinate range. - resolve_support / select_seeds: per-label k-mer support of a query (explicit labels or discovery by column or header), seed candidates, and the max_support / longest_first seed selection. - CoordToHeader::find_header (a lazily built reverse header index owned by the object) and last_coord. Source paths on d5d26baa: metagraph/src/graph/traversal/{traversal_types,label_oracle,resolve}.{hpp,cpp} metagraph/src/annotation/coord_to_header.{hpp,cpp} metagraph/tests/graph/traversal/{test_label_oracle,test_resolve}.cpp Co-authored-by: Gunnar Rätsch Co-Authored-By: Claude Opus 5.5 --- metagraph/src/annotation/coord_to_header.cpp | 50 ++ metagraph/src/annotation/coord_to_header.hpp | 27 +- .../src/graph/traversal/label_oracle.cpp | 280 ++++++++++ .../src/graph/traversal/label_oracle.hpp | 167 ++++++ metagraph/src/graph/traversal/resolve.cpp | 525 ++++++++++++++++++ metagraph/src/graph/traversal/resolve.hpp | 160 ++++++ .../src/graph/traversal/traversal_types.cpp | 62 +++ .../src/graph/traversal/traversal_types.hpp | 74 +++ .../graph/traversal/test_label_oracle.cpp | 377 +++++++++++++ .../tests/graph/traversal/test_resolve.cpp | 403 ++++++++++++++ 10 files changed, 2124 insertions(+), 1 deletion(-) create mode 100644 metagraph/src/graph/traversal/label_oracle.cpp create mode 100644 metagraph/src/graph/traversal/label_oracle.hpp create mode 100644 metagraph/src/graph/traversal/resolve.cpp create mode 100644 metagraph/src/graph/traversal/resolve.hpp create mode 100644 metagraph/src/graph/traversal/traversal_types.cpp create mode 100644 metagraph/src/graph/traversal/traversal_types.hpp create mode 100644 metagraph/tests/graph/traversal/test_label_oracle.cpp create mode 100644 metagraph/tests/graph/traversal/test_resolve.cpp diff --git a/metagraph/src/annotation/coord_to_header.cpp b/metagraph/src/annotation/coord_to_header.cpp index 272a374f54..e03ddec80f 100644 --- a/metagraph/src/annotation/coord_to_header.cpp +++ b/metagraph/src/annotation/coord_to_header.cpp @@ -1,5 +1,6 @@ #include "annotation/coord_to_header.hpp" +#include #include #include "common/logger.hpp" @@ -35,6 +36,37 @@ struct HeaderOfRecordEqual { } // namespace +// header -> (column, seq_id) of its first occurrence in column order +struct CoordToHeader::HeaderIndex { + tsl::hopscotch_map> map; +}; + +CoordToHeader::CoordToHeader() {} +CoordToHeader::~CoordToHeader() {} + +std::optional> +CoordToHeader::find_header(std::string_view header) const { + std::lock_guard lock(header_index_mutex_); + if (!header_index_) { + auto index = std::make_unique(); + size_t total = 0; + for (const auto &column : headers_) { + total += column.size(); + } + index->map.reserve(total); + for (Column col = 0; col < headers_.size(); ++col) { + for (size_t s = 0; s < headers_[col].size(); ++s) { + index->map.try_emplace(std::string_view(headers_[col][s]), col, s); + } + } + header_index_ = std::move(index); + } + auto it = header_index_->map.find(header); + if (it == header_index_->map.end()) + return std::nullopt; + return it->second; +} + void CoordToHeader::find_duplicate_headers(const std::vector &headers, const std::function &callback) { tsl::hopscotch_set @@ -103,6 +135,19 @@ CoordToHeader::map_single_coord(Column col, uint64_t coord) const { return { header, local_coord }; } +uint64_t CoordToHeader::last_coord(Column col, size_t seq_id) const { + if (col >= num_columns()) { + throw std::out_of_range(fmt::format("Column {} out of range " + "(CoordToHeader has {} columns)", col, num_columns())); + } + if (seq_id >= num_sequences(col)) { + throw std::out_of_range(fmt::format("Sequence id {} out of range for column {} " + "({} sequences)", seq_id, col, num_sequences(col))); + } + // a set bit marks the last coordinate of each sequence + return coord_offsets_[col].select1(seq_id + 1); +} + uint64_t CoordToHeader::num_kmers_in_sequence(Column col, size_t seq_id) const { if (col >= num_columns()) { throw std::out_of_range(fmt::format("Column {} out of range " @@ -129,6 +174,11 @@ bool CoordToHeader::load(const std::string &filename_base) { utils::file_read_failure_detail(path)); return false; } + { + // the index views the headers about to be replaced + std::lock_guard lock(header_index_mutex_); + header_index_.reset(); + } try { uint64_t num_columns = load_number(*in); diff --git a/metagraph/src/annotation/coord_to_header.hpp b/metagraph/src/annotation/coord_to_header.hpp index 34c8a61d09..6fe0f72651 100644 --- a/metagraph/src/annotation/coord_to_header.hpp +++ b/metagraph/src/annotation/coord_to_header.hpp @@ -2,7 +2,11 @@ #define __COORD_TO_HEADER_HPP__ #include +#include +#include +#include #include +#include #include #include "annotation/int_matrix/base/int_matrix.hpp" @@ -42,10 +46,12 @@ class CoordToHeader { using Tuple = mtg::annot::matrix::MultiIntMatrix::Tuple; using RowTuples = mtg::annot::matrix::MultiIntMatrix::RowTuples; - CoordToHeader() {} + // out of line, as the destructor: the header index type is private to the .cpp + CoordToHeader(); // Constructor from sequence names and number of k-mers in them, per annotation column. CoordToHeader(std::vector> &&headers, std::vector> &&num_kmers); + ~CoordToHeader(); bool load(const std::string &filename_base); void serialize(const std::string &filename_base) const; @@ -53,6 +59,10 @@ class CoordToHeader { // global_coord -> (seq_id, local_coord) std::pair map_single_coord(Column col, uint64_t coord) const; + // The last global coordinate of a sequence (its first is the previous sequence's + // last + 1). Throws as num_kmers_in_sequence does. + uint64_t last_coord(Column col, size_t seq_id) const; + /** * Batch variant of `map_single_coord` that replaces each global coord * in `rows_tuples` with a packed `(seq_id, local_coord)` pair @@ -78,6 +88,14 @@ class CoordToHeader { // FASTA headers of the sequences in a column, indexed by seq_id. const std::vector& get_headers(Column column) const { return headers_[column]; } + /** + * The (column, seq_id) of a FASTA header: its first occurrence in column order (a + * header is unique within a column, not within the index), or nothing. The reverse + * index is built on the first call and owned by this object, so it lives and dies with + * the headers it indexes. load() drops it. Thread-safe. + */ + std::optional> find_header(std::string_view header) const; + /** * The records of |headers| whose header repeats an earlier record's, in record order: * |callback|(first, repeat) with the 0-based indices of the first occurrence and of the @@ -91,8 +109,15 @@ class CoordToHeader { static constexpr auto kExtension = ".seqs"; private: + struct HeaderIndex; + std::vector> headers_; std::vector coord_offsets_; + + // find_header(): keys are views into |headers_|, which only load() changes, and load() + // drops the index first + mutable std::mutex header_index_mutex_; + mutable std::unique_ptr header_index_; }; } // namespace annot diff --git a/metagraph/src/graph/traversal/label_oracle.cpp b/metagraph/src/graph/traversal/label_oracle.cpp new file mode 100644 index 0000000000..60f1c89236 --- /dev/null +++ b/metagraph/src/graph/traversal/label_oracle.cpp @@ -0,0 +1,280 @@ +#include "label_oracle.hpp" + +#include +#include + +#include "graph/annotated_dbg.hpp" +#include "graph/representation/canonical_dbg.hpp" +#include "graph/representation/succinct/dbg_succinct.hpp" +#include "graph/graph_extensions/node_first_cache.hpp" +#include "annotation/coord_to_header.hpp" +#include "common/logger.hpp" +#include "common/unix_tools.hpp" + + +namespace mtg { +namespace graph { +namespace traversal { + +using mtg::common::logger; +using annot::matrix::BinaryMatrix; +using annot::matrix::MultiIntMatrix; + + +LabelOracle::LabelOracle(const AnnotatedDBG &anno_graph, + const annot::CoordToHeader *coord_to_header_override) + : anno_graph_(anno_graph) { + const DeBruijnGraph &loaded = anno_graph_.get_graph(); + + // A CanonicalDBG wrapper means a PRIMARY index. Otherwise the mode tells native + // CANONICAL graphs from BASIC ones. + if (const auto *wrapper = dynamic_cast(&loaded)) { + regime_ = Regime::PRIMARY; + // a wrapper of our own: same node ids, own caches; the shared one is never mutated + local_canonical_ = std::make_shared(wrapper->get_graph_ptr()); + if (const auto *dbg_succ = dynamic_cast(&wrapper->get_graph())) + local_canonical_->add_extension(std::make_shared(*dbg_succ)); + graph_ = local_canonical_.get(); + } else { + if (loaded.get_mode() == DeBruijnGraph::PRIMARY) { + throw std::invalid_argument("PRIMARY graphs must be wrapped into CanonicalDBG " + "before traversal"); + } + regime_ = loaded.get_mode() == DeBruijnGraph::CANONICAL ? Regime::CANONICAL + : Regime::BASIC; + graph_ = &loaded; + } + key_mapper_ = std::make_unique(*graph_); + + // get_matrix() may lock and flush (ColumnCompressed): call it once per request + matrix_ = &anno_graph_.get_annotator().get_matrix(); + tuples_ = dynamic_cast(matrix_); + coord_to_header_ = coord_to_header_override ? coord_to_header_override + : anno_graph_.get_coord_to_header(); + if (coord_to_header_ && !tuples_) { + logger->warn("CoordToHeader present but the annotation has no coordinates; " + "header labels are unavailable"); + coord_to_header_ = nullptr; + } +} + +LabelOracle::~LabelOracle() {} + +size_t LabelOracle::get_k() const { return graph_->get_k(); } + +uint64_t LabelOracle::num_columns() const { return matrix_->num_columns(); } + +std::vector LabelOracle::keys_of_sequence(std::string_view sequence) const { + std::vector keys; + if (sequence.size() < get_k()) + return keys; + keys.reserve(sequence.size() - get_k() + 1); + // map_to_nodes returns annotation keys in every regime: the node itself (BASIC), + // the primary node (CanonicalDBG), the canonical node (native CANONICAL graphs) + graph_->map_to_nodes(sequence, [&](node_index key) { keys.push_back(key); }); + counters_.keys_mapped += keys.size(); + return keys; +} + +node_index LabelOracle::key_of(node_index node, std::string_view kmer) const { + counters_.keys_mapped++; + return node == npos ? npos : key_mapper_->key(node, kmer); +} + +std::optional LabelOracle::find_column(const std::string &name) const { + const auto &encoder = anno_graph_.get_annotator().get_label_encoder(); + if (!encoder.label_exists(name)) + return std::nullopt; + return static_cast(encoder.encode(name)); +} + +const std::string& LabelOracle::column_name(Column column) const { + return anno_graph_.get_annotator().get_label_encoder().decode(column); +} + +std::optional LabelOracle::find_header(const std::string &name) const { + if (!coord_to_header_) + return std::nullopt; + auto found = coord_to_header_->find_header(name); + if (!found) + return std::nullopt; + LabelRef ref; + ref.kind = LabelKind::HEADER; + ref.column = found->first; + ref.seq_id = found->second; + ref.name = name; + return ref; +} + +const std::string& LabelOracle::header_name(Column column, uint64_t seq_id) const { + assert(coord_to_header_); + return coord_to_header_->get_headers(column).at(seq_id); +} + +std::pair LabelOracle::map_coord(Column column, uint64_t coord) const { + assert(coord_to_header_); + counters_.coords_mapped++; + return coord_to_header_->map_single_coord(column, coord); +} + +std::pair LabelOracle::sequence_coords(Column column, uint64_t seq_id) const { + assert(coord_to_header_); + const uint64_t last = coord_to_header_->last_coord(column, seq_id); + const uint64_t first = seq_id ? coord_to_header_->last_coord(column, seq_id - 1) + 1 : 0; + return { first, last }; +} + +LabelRef LabelOracle::resolve_label(const std::string &name) const { + if (auto column = find_column(name)) { + LabelRef ref; + ref.kind = LabelKind::COLUMN; + ref.column = *column; + ref.name = name; + return ref; + } + if (auto header = find_header(name)) + return *header; + throw std::invalid_argument("Label not found in the annotation: '" + name + "'" + + (coord_to_header_ ? "" : " (no sequence header index loaded)")); +} + +std::vector +LabelOracle::get_rows(const std::vector &rows) const { + Timer timer; + auto result = matrix_->get_rows(rows); + counters_.rows_fetched += rows.size(); + counters_.fetch_seconds += timer.elapsed(); + return result; +} + +std::vector +LabelOracle::get_row_tuples(const std::vector &rows) const { + if (!tuples_) + throw std::logic_error("The annotation has no coordinates"); + Timer timer; + auto result = tuples_->get_row_tuples(rows); + counters_.tuple_rows_fetched += rows.size(); + counters_.fetch_seconds += timer.elapsed(); + return result; +} + + +LabelQuery::LabelQuery(const LabelOracle &oracle, + std::vector labels, + size_t max_cache_size) + : oracle_(oracle), labels_(std::move(labels)), max_cache_size_(max_cache_size) { + if (labels_.empty()) + throw std::invalid_argument("Empty permitted label set"); + + std::set> headers; + for (LabelId i = 0; i < labels_.size(); ++i) { + const LabelRef &ref = labels_[i]; + if (ref.kind == LabelKind::COLUMN) { + if (!column_labels_.try_emplace(ref.column, i).second) + throw std::invalid_argument("Duplicate label '" + ref.name + "'"); + continue; + } + if (!oracle_.coord_to_header()) + throw std::invalid_argument("Sequence header labels require an annotation with " + "coordinates and a CoordToHeader index"); + if (!headers.emplace(ref.column, ref.seq_id).second) + throw std::invalid_argument("Duplicate label '" + ref.name + "'"); + const auto [first, last] = oracle_.sequence_coords(ref.column, ref.seq_id); + header_ranges_[ref.column].push_back(HeaderRange{ first, last, i }); + } + for (auto it = header_ranges_.begin(); it != header_ranges_.end(); ++it) { + std::sort(it.value().begin(), it.value().end(), + [](const HeaderRange &a, const HeaderRange &b) { return a.first < b.first; }); + } + tuples_ = !header_ranges_.empty(); +} + +void LabelQuery::hits_from_row(const BinaryMatrix::SetBitPositions &row, NodeHits *hits) const { + for (Column c : row) { + auto it = column_labels_.find(c); + if (it != column_labels_.end()) + hits->push_back(it->second); + } + std::sort(hits->begin(), hits->end()); +} + +void LabelQuery::hits_from_tuples(const MultiIntMatrix::RowTuples &row, NodeHits *hits) const { + for (const auto &[c, coords] : row) { + if (auto it = column_labels_.find(c); it != column_labels_.end()) + hits->push_back(it->second); + auto rt = header_ranges_.find(c); + if (rt == header_ranges_.end()) + continue; + const std::vector &ranges = rt->second; + for (uint64_t coord : coords) { + // the last range starting at or before the coordinate (the ranges are disjoint) + auto it = std::upper_bound(ranges.begin(), ranges.end(), coord, + [](uint64_t x, const HeaderRange &r) { return x < r.first; }); + if (it != ranges.begin() && coord <= (--it)->last) + hits->push_back(it->label); + } + } + std::sort(hits->begin(), hits->end()); + hits->erase(std::unique(hits->begin(), hits->end()), hits->end()); +} + +void LabelQuery::fetch_uncached(const std::vector &keys) { + if (keys.empty()) + return; + if (cache_.size() + keys.size() > max_cache_size_) + cache_.clear(); + + std::vector rows; + rows.reserve(keys.size()); + for (node_index key : keys) { + rows.push_back(AnnotatedDBG::graph_to_anno_index(key)); + } + if (tuples_) { + auto tuple_rows = oracle_.get_row_tuples(rows); + for (size_t i = 0; i < keys.size(); ++i) { + hits_from_tuples(tuple_rows[i], &cache_[keys[i]]); + } + } else { + auto set_rows = oracle_.get_rows(rows); + for (size_t i = 0; i < keys.size(); ++i) { + hits_from_row(set_rows[i], &cache_[keys[i]]); + } + } +} + +std::vector LabelQuery::fetch(const std::vector &keys) { + auto &counters = oracle_.counters(); + std::vector missing; + tsl::hopscotch_set requested; + for (node_index key : keys) { + if (key == npos) + continue; + counters.rows_requested++; + if (cache_.count(key)) { + counters.cache_hits++; + } else if (requested.insert(key).second) { + missing.push_back(key); + } + } + // a clear in fetch_uncached may drop keys requested here: read them all again then + if (cache_.size() + missing.size() > max_cache_size_) { + missing.clear(); + requested.clear(); + for (node_index key : keys) { + if (key != npos && requested.insert(key).second) + missing.push_back(key); + } + } + fetch_uncached(missing); + + std::vector result(keys.size()); + for (size_t i = 0; i < keys.size(); ++i) { + if (keys[i] != npos) + result[i] = cache_.at(keys[i]); + } + return result; +} + +} // namespace traversal +} // namespace graph +} // namespace mtg diff --git a/metagraph/src/graph/traversal/label_oracle.hpp b/metagraph/src/graph/traversal/label_oracle.hpp new file mode 100644 index 0000000000..a92cf5e224 --- /dev/null +++ b/metagraph/src/graph/traversal/label_oracle.hpp @@ -0,0 +1,167 @@ +#ifndef __TRAVERSAL_LABEL_ORACLE_HPP__ +#define __TRAVERSAL_LABEL_ORACLE_HPP__ + +#include +#include +#include +#include +#include + +#include +#include + +#include "traversal_types.hpp" +#include "graph/annotation_key_mapper.hpp" +#include "annotation/int_matrix/base/int_matrix.hpp" + + +namespace mtg { + +namespace annot { +class CoordToHeader; +} // namespace annot + +namespace graph { + +class AnnotatedDBG; +class DeBruijnGraph; +class CanonicalDBG; + +namespace traversal { + +/** + * Per-request view of an annotated graph for traversal: the oriented graph to walk, + * the mapping of its nodes to annotation rows, label names, and batched row access. + * + * In the PRIMARY regime the oracle walks a CanonicalDBG of its own over the loaded primary + * graph (same node ids, its own caches), so concurrent requests never share wrapper state. + * Not thread-safe: one oracle per request. + */ +class LabelOracle { + public: + struct Counters { + uint64_t keys_mapped = 0; // oriented nodes mapped to annotation keys + uint64_t rows_requested = 0; // annotation rows requested (incl. cache hits) + uint64_t cache_hits = 0; + uint64_t rows_fetched = 0; // rows reconstructed with get_rows + uint64_t tuple_rows_fetched = 0; // rows reconstructed with get_row_tuples + uint64_t coords_mapped = 0; // coordinates mapped to (sequence, local coord) + double fetch_seconds = 0; + }; + + // |coord_to_header_override| replaces the mapping loaded with the index (tests). + explicit LabelOracle(const AnnotatedDBG &anno_graph, + const annot::CoordToHeader *coord_to_header_override = nullptr); + ~LabelOracle(); + + const AnnotatedDBG& anno_graph() const { return anno_graph_; } + // the graph used for traversal (oriented node ids) + const DeBruijnGraph& graph() const { return *graph_; } + Regime regime() const { return regime_; } + size_t get_k() const; + + bool has_coordinates() const { return tuples_; } + // the CoordToHeader index, or nullptr when there is none or the annotation has no + // coordinates (header labels then do not exist) + const annot::CoordToHeader* coord_to_header() const { return coord_to_header_; } + uint64_t num_columns() const; + + // Annotation keys of all k-mers of |sequence| (npos for k-mers not in the graph). + std::vector keys_of_sequence(std::string_view sequence) const; + // Annotation key of the oriented |node| spelling |kmer|. + node_index key_of(node_index node, std::string_view kmer) const; + + // Label resolution + std::optional find_column(const std::string &name) const; + const std::string& column_name(Column column) const; + // The first (column, sequence) holding a FASTA header (CoordToHeader::find_header) + std::optional find_header(const std::string &name) const; + const std::string& header_name(Column column, uint64_t seq_id) const; + // (seq_id, local coordinate) of a column coordinate. Requires a CoordToHeader. + std::pair map_coord(Column column, uint64_t coord) const; + // The first and last column coordinates of sequence |seq_id|. Requires a CoordToHeader. + std::pair sequence_coords(Column column, uint64_t seq_id) const; + // A column name first, then a header. Throws std::invalid_argument naming the label if + // neither exists. + LabelRef resolve_label(const std::string &name) const; + + // Batched row access (rows = annotation key - 1). Updates the counters. + std::vector + get_rows(const std::vector &rows) const; + std::vector + get_row_tuples(const std::vector &rows) const; + + Counters& counters() const { return counters_; } + + private: + const AnnotatedDBG &anno_graph_; + std::shared_ptr local_canonical_; + const DeBruijnGraph *graph_ = nullptr; + std::unique_ptr key_mapper_; + Regime regime_ = Regime::BASIC; + + const annot::matrix::BinaryMatrix *matrix_ = nullptr; + const annot::matrix::MultiIntMatrix *tuples_ = nullptr; + const annot::CoordToHeader *coord_to_header_ = nullptr; + + mutable Counters counters_; +}; + + +/** + * Batched membership of a fixed set of permitted labels at annotation keys. Column + * labels are read from whole rows; header labels (one sequence of a column) from tuple + * rows, a coordinate counting for a label when it falls into that label's sequence. + * Results are cached per key, up to |max_cache_size| keys, then the cache is cleared. + */ +class LabelQuery { + public: + // the permitted labels present at a key, ascending + using NodeHits = std::vector; + + // Throws std::invalid_argument for an empty or duplicate label set, or header labels + // on an annotation without coordinates. + LabelQuery(const LabelOracle &oracle, + std::vector labels, + size_t max_cache_size = 1'000'000); + + const std::vector& labels() const { return labels_; } + // which accessor is used: "rows" or "tuples" + const char* access_path() const { return tuples_ ? "tuples" : "rows"; } + + // Hits for every annotation key (npos yields no hits). + std::vector fetch(const std::vector &keys); + + private: + struct HeaderRange { + uint64_t first; + uint64_t last; + LabelId label; + }; + + void fetch_uncached(const std::vector &keys); + void hits_from_row(const annot::matrix::BinaryMatrix::SetBitPositions &row, + NodeHits *hits) const; + void hits_from_tuples(const annot::matrix::MultiIntMatrix::RowTuples &row, + NodeHits *hits) const; + + const LabelOracle &oracle_; + std::vector labels_; + bool tuples_ = false; + size_t max_cache_size_; + + // column -> label id for COLUMN labels + tsl::hopscotch_map column_labels_; + // per column holding HEADER labels, the coordinate ranges of the requested sequences, + // ascending: a coordinate is looked up among them by a binary search instead of being + // mapped to its sequence + tsl::hopscotch_map> header_ranges_; + + tsl::hopscotch_map cache_; +}; + +} // namespace traversal +} // namespace graph +} // namespace mtg + +#endif // __TRAVERSAL_LABEL_ORACLE_HPP__ diff --git a/metagraph/src/graph/traversal/resolve.cpp b/metagraph/src/graph/traversal/resolve.cpp new file mode 100644 index 0000000000..e1bd6bf25f --- /dev/null +++ b/metagraph/src/graph/traversal/resolve.cpp @@ -0,0 +1,525 @@ +#include "resolve.hpp" + +#include +#include +#include + +#include + +#include "graph/annotated_dbg.hpp" +#include "annotation/coord_to_header.hpp" +#include "common/seq_tools/reverse_complement.hpp" + + +namespace mtg { +namespace graph { +namespace traversal { + +using annot::matrix::BinaryMatrix; +using annot::matrix::MultiIntMatrix; + +uint64_t fnv1a64(std::string_view data, uint64_t hash) { + for (unsigned char c : data) { + hash ^= c; + hash *= 0x100000001b3ULL; + } + return hash; +} + +std::string hex64(uint64_t x) { + static const char *digits = "0123456789abcdef"; + std::string s(16, '0'); + for (int i = 15; i >= 0; --i) { + s[i] = digits[x & 15]; + x >>= 4; + } + return s; +} + +std::string make_seed_id(std::string_view sequence, + bool canonical_orientation, + std::vector labels) { + // defined over the sequence after the build's case mapping, so that a seed resubmitted + // in another case keeps its id + std::string seq(sequence); +#if ! _DNA_CASE_SENSITIVE_GRAPH + for (char &c : seq) { + c = std::toupper(static_cast(c)); + } +#endif + if (canonical_orientation) { + std::string rc(seq); + ::reverse_complement(rc); + if (rc < seq) + seq = rc; + } + std::sort(labels.begin(), labels.end()); + // an empty release name and a tab come first: the id of the labeled-traversal + // prototype for an index without a release name + uint64_t h = fnv1a64("\t"); + h = fnv1a64(seq, h); + h = fnv1a64("\t", h); + for (const auto &l : labels) { + h = fnv1a64(std::to_string(l.size()) + ":" + l, h); + } + return hex64(h); +} + +std::string encode_runs(const std::vector &runs, uint64_t num_kmers) { + std::string out; + uint64_t pos = 0; + for (const auto &run : runs) { + if (run.begin > pos) + out += "o" + std::to_string(run.begin - pos); + out += "x" + std::to_string(run.size()); + pos = run.end; + } + if (num_kmers > pos) + out += "o" + std::to_string(num_kmers - pos); + return out; +} + + +namespace { + +// The support of one label, built k-mer by k-mer in ascending order (each k-mer at most once) +struct LabelSupport { + uint64_t kmers = 0; + std::vector runs; + + void add(uint64_t i) { + // present at the k-mer once, however many hits name it + if (!runs.empty() && runs.back().end == i + 1) + return; + kmers++; + if (!runs.empty() && runs.back().end == i) { + runs.back().end = i + 1; + } else { + runs.push_back({ i, i + 1 }); + } + } +}; + +std::vector anno_rows(const std::vector &keys) { + std::vector rows; + rows.reserve(keys.size()); + for (node_index key : keys) { + rows.push_back(AnnotatedDBG::graph_to_anno_index(key)); + } + return rows; +} + +// The labels of |options.discover_kind| present at the k-mers |pos| (keys |keys|), each with +// its support, in the order first seen. Rows are decoded kResolveBatchRows at a time. +void discover_labels(const LabelOracle &oracle, + const std::vector &keys, + const std::vector &pos, + LabelKind kind, + std::vector> *ids, + std::vector *found) { + std::map, size_t> index; + auto of = [&](Column c, uint64_t seq_id) -> LabelSupport& { + auto [it, inserted] = index.try_emplace({ c, seq_id }, found->size()); + if (inserted) { + found->emplace_back(); + ids->emplace_back(c, seq_id); + } + return (*found)[it->second]; + }; + for (size_t begin = 0; begin < keys.size(); begin += kResolveBatchRows) { + const size_t end = std::min(keys.size(), begin + kResolveBatchRows); + const std::vector batch(keys.begin() + begin, keys.begin() + end); + if (kind == LabelKind::COLUMN) { + const auto rows = oracle.get_rows(anno_rows(batch)); + for (size_t j = 0; j < rows.size(); ++j) { + for (Column c : rows[j]) { + of(c, 0).add(pos[begin + j]); + } + } + } else { + const auto rows = oracle.get_row_tuples(anno_rows(batch)); + for (size_t j = 0; j < rows.size(); ++j) { + for (const auto &[c, coords] : rows[j]) { + for (uint64_t coord : coords) { + of(c, oracle.map_coord(c, coord).first).add(pos[begin + j]); + } + } + } + } + } +} + +} // namespace + +SupportProfile resolve_support(const LabelOracle &oracle, + std::string_view query, + const ResolveOptions &options) { + SupportProfile profile; + profile.k = oracle.get_k(); + profile.regime = oracle.regime(); + + if (options.labels.empty() == !options.discover) + throw std::invalid_argument("Specify either an explicit label list or discovery, not both"); + if (query.size() < profile.k) + throw std::invalid_argument("Query shorter than k"); + + profile.num_kmers = query.size() - profile.k + 1; + const std::vector keys = oracle.keys_of_sequence(query); + assert(keys.size() == profile.num_kmers); + + std::vector present_keys; + std::vector present_pos; // the k-mer index of each present key + for (uint64_t i = 0; i < keys.size(); ++i) { + if (keys[i] != npos) { + present_keys.push_back(keys[i]); + present_pos.push_back(i); + } + } + profile.graph_runs = runs_of(profile.num_kmers, [&](uint64_t i) { return keys[i] != npos; }); + + std::vector refs; + std::vector support; + if (!options.labels.empty()) { + for (const auto &name : options.labels) { + refs.push_back(oracle.resolve_label(name)); + } + LabelQuery query_labels(oracle, refs); + const auto hits = query_labels.fetch(keys); + support.resize(refs.size()); + for (uint64_t i = 0; i < hits.size(); ++i) { + for (LabelId l : hits[i]) { + support[l].add(i); + } + } + } else { + if (options.discover_kind == LabelKind::HEADER && !oracle.coord_to_header()) + throw std::invalid_argument("Header discovery requires a CoordToHeader index"); + std::vector> ids; + std::vector found; + discover_labels(oracle, present_keys, present_pos, options.discover_kind, &ids, &found); + + // more k-mers first; ties by column id, then seq_id + std::vector ranked(found.size()); + for (size_t i = 0; i < ranked.size(); ++i) { + ranked[i] = i; + } + std::sort(ranked.begin(), ranked.end(), [&](size_t a, size_t b) { + return found[a].kmers != found[b].kmers ? found[a].kmers > found[b].kmers + : ids[a] < ids[b]; + }); + if (ranked.size() > options.discover_max_labels) { + LabelTruncation trunc; + trunc.total = ranked.size(); + trunc.kept = options.discover_max_labels; + trunc.min_kept_kmers = found[ranked[trunc.kept - 1]].kmers; + trunc.max_dropped_kmers = found[ranked[trunc.kept]].kmers; + for (size_t i = trunc.kept; i < ranked.size(); ++i) { + if (found[ranked[i]].kmers == present_keys.size()) + trunc.dropped_full_length++; + } + ranked.resize(trunc.kept); + profile.labels_truncated = trunc; + } + for (size_t at : ranked) { + LabelRef ref; + ref.kind = options.discover_kind; + ref.column = ids[at].first; + ref.seq_id = ids[at].second; + ref.name = ref.kind == LabelKind::COLUMN ? oracle.column_name(ref.column) + : oracle.header_name(ref.column, ref.seq_id); + refs.push_back(std::move(ref)); + support.push_back(std::move(found[at])); + } + } + + for (LabelId l = 0; l < refs.size(); ++l) { + LabelProfile lp; + lp.label = refs[l]; + lp.kmers_supported = support[l].kmers; + lp.runs = std::move(support[l].runs); + profile.labels.push_back(std::move(lp)); + } + + // ---- seed candidates: identical maximal runs grouped + std::map> groups; + for (LabelId l = 0; l < profile.labels.size(); ++l) { + for (const auto &run : profile.labels[l].runs) { + if (run.size() >= options.min_block_kmers) + groups[run].push_back(l); + } + } + for (auto &[iv, labels] : groups) { + profile.candidates.push_back({ iv, std::move(labels) }); + } + std::sort(profile.candidates.begin(), profile.candidates.end(), + [](const SeedCandidate &a, const SeedCandidate &b) { + if (a.kmers.size() != b.kmers.size()) + return a.kmers.size() > b.kmers.size(); + if (a.labels.size() != b.labels.size()) + return a.labels.size() > b.labels.size(); + return a.kmers.begin < b.kmers.begin; + }); + return profile; +} + + +namespace { + +// A Fenwick tree of counts over positions 0..n-1: add, the sum over positions < i, and +// the position of the k-th unit, each in O(log n). +struct Fenwick { + std::vector tree; + size_t n; + size_t top = 1; // the largest power of two <= n (1 when n == 0) + + explicit Fenwick(size_t size) : tree(size + 1, 0), n(size) { + while (top * 2 <= n) { + top *= 2; + } + } + void add(size_t i, uint32_t v) { + for (++i; i <= n; i += i & (~i + 1)) { + tree[i] += v; + } + } + uint32_t prefix(size_t i) const { // sum over positions < i + uint32_t s = 0; + for (; i > 0; i -= i & (~i + 1)) { + s += tree[i]; + } + return s; + } + size_t kth(uint32_t k) const { // the 0-based position holding the k-th unit (1-based k) + size_t pos = 0; + for (size_t step = top; step; step >>= 1) { + if (pos + step <= n && tree[pos + step] < k) { + pos += step; + k -= tree[pos]; + } + } + return pos; + } +}; + +// does label l have a run containing [begin, end)? +bool covers(const LabelProfile &lp, const KmerInterval &iv) { + for (const auto &run : lp.runs) { + if (run.begin <= iv.begin && run.end >= iv.end) + return true; + } + return false; +} + +std::vector covering_labels(const SupportProfile &profile, const KmerInterval &iv) { + std::vector out; + for (LabelId l = 0; l < profile.labels.size(); ++l) { + if (covers(profile.labels[l], iv)) + out.push_back(l); + } + return out; +} + +struct Picked { + KmerInterval kmers; + std::vector labels; +}; + +// MAX_SUPPORT: repeatedly pick the interval of at least |min_kmers| k-mers contained in the +// most support runs, outside the intervals picked before. +// +// For a fixed start a, the count |{runs [s, e) : s <= a, e >= b}| is non-increasing in b, so +// the best count is reached at the shortest admissible b = a + min_kmers, and the longest +// interval with that count ends at b = min{e : e >= a + min_kmers}. So one search per +// candidate a suffices, and a ranges over run starts and over (e - min_kmers): O(R log R) per +// round for R runs. +std::vector pick_max_support(const SupportProfile &profile, uint64_t min_kmers, + size_t max_seeds) { + struct Run { uint64_t begin, end; }; + std::vector all_runs; + for (const auto &lp : profile.labels) { + for (const auto &run : lp.runs) { + if (run.size() >= min_kmers) + all_runs.push_back({ run.begin, run.end }); + } + } + std::sort(all_runs.begin(), all_runs.end(), + [](const Run &a, const Run &b) { return a.begin < b.begin; }); + + std::vector picked; + // intervals already claimed by a picked seed, as sorted disjoint ranges + std::vector taken; + std::vector runs; + while (picked.size() < max_seeds) { + // the runs with the claimed intervals subtracted, so that later rounds see the + // clipped runs' own endpoints (with A=[0,10), B=[0,5) and [0,5) taken, [5,10) is + // only reachable through them) + runs.clear(); + for (const Run &run : all_runs) { + uint64_t begin = run.begin; + for (const auto &t : taken) { + if (t.end <= begin || t.begin >= run.end) + continue; + if (t.begin > begin && t.begin - begin >= min_kmers) + runs.push_back({ begin, t.begin }); + begin = std::max(begin, t.end); + } + if (begin < run.end && run.end - begin >= min_kmers) + runs.push_back({ begin, run.end }); + } + std::sort(runs.begin(), runs.end(), + [](const Run &a, const Run &b) { return a.begin < b.begin; }); + if (runs.empty()) + break; + + std::vector candidate_a; + for (const auto &run : runs) { + candidate_a.push_back(run.begin); + if (run.end >= min_kmers) + candidate_a.push_back(run.end - min_kmers); + } + std::sort(candidate_a.begin(), candidate_a.end()); + candidate_a.erase(std::unique(candidate_a.begin(), candidate_a.end()), candidate_a.end()); + + // the ends of the runs with begin <= a, as counts over the compressed run ends + std::vector ends_at; + for (const Run &run : runs) { + ends_at.push_back(run.end); + } + std::sort(ends_at.begin(), ends_at.end()); + ends_at.erase(std::unique(ends_at.begin(), ends_at.end()), ends_at.end()); + auto index_of = [&](uint64_t x) { + return static_cast(std::lower_bound(ends_at.begin(), ends_at.end(), x) + - ends_at.begin()); + }; + Fenwick ends(ends_at.size()); + uint32_t inserted = 0; + size_t next_run = 0; + std::optional best; + size_t best_count = 0; + for (uint64_t a : candidate_a) { + while (next_run < runs.size() && runs[next_run].begin <= a) { + ends.add(index_of(runs[next_run].end), 1); + inserted++; + next_run++; + } + if (a + min_kmers > profile.num_kmers) + break; + const size_t lo = index_of(a + min_kmers); + const uint32_t below = lo < ends_at.size() ? ends.prefix(lo) : inserted; + if (below == inserted) + continue; + const KmerInterval iv { a, ends_at[ends.kth(below + 1)] }; + const size_t count = inserted - below; + if (!best || count > best_count + || (count == best_count && (iv.size() > best->size() + || (iv.size() == best->size() && iv.begin < best->begin)))) { + best = iv; + best_count = count; + } + } + if (!best) + break; + picked.push_back({ *best, covering_labels(profile, *best) }); + assert(picked.back().labels.size() == best_count); + taken.push_back(*best); + std::sort(taken.begin(), taken.end()); + } + return picked; +} + +} // namespace + +SeedSelection select_seeds(const SupportProfile &profile, + std::string_view query, + const SelectionPolicy &policy_in, + bool canonical_orientation) { + SeedSelection selection; + selection.policy = policy_in; + SelectionPolicy &policy = selection.policy; + if (policy.min_block_bp < profile.k) + policy.min_block_bp = profile.k; + const uint64_t min_kmers = policy.min_block_bp - profile.k + 1; + selection.num_candidates = profile.candidates.size(); + + std::vector picked; + if (policy.policy == SelectionPolicy::LONGEST_FIRST) { + for (const auto &cand : profile.candidates) { + if (cand.kmers.size() < min_kmers) + continue; + selection.num_eligible++; + if (picked.size() < policy.max_seeds) + picked.push_back({ cand.kmers, cand.labels }); + } + } else { + picked = pick_max_support(profile, min_kmers, policy.max_seeds); + selection.num_eligible = picked.size(); + } + + // merge overlapping picks onto their common interval + if (policy.merge_overlapping) { + std::sort(picked.begin(), picked.end(), + [](const Picked &a, const Picked &b) { return a.kmers < b.kmers; }); + std::vector merged; + for (auto &p : picked) { + if (!merged.empty() && merged.back().kmers.end > p.kmers.begin) { + Picked &m = merged.back(); + KmerInterval common { std::max(m.kmers.begin, p.kmers.begin), + std::min(m.kmers.end, p.kmers.end) }; + if (common.size() >= min_kmers) { + m.kmers = common; + m.labels = covering_labels(profile, common); + continue; + } + } + merged.push_back(std::move(p)); + } + picked.swap(merged); + } + + // freeze + for (auto &p : picked) { + FrozenSeed seed; + seed.kmers = p.kmers; + seed.sequence = std::string(query.substr(p.kmers.begin, p.kmers.size() + profile.k - 1)); + + // the labels in the order of the hash of their name (then the name) + std::vector> ordered; + for (LabelId l : p.labels) { + const std::string &name = profile.labels[l].label.name; + ordered.emplace_back(fnv1a64(name + "\t0"), name); + } + std::sort(ordered.begin(), ordered.end()); + + seed.population.supporting_total = ordered.size(); + uint64_t digest = 0xcbf29ce484222325ULL; + for (size_t i = 0; i < ordered.size(); ++i) { + if (i < policy.max_labels_per_seed) { + seed.labels.push_back(ordered[i].second); + } else { + seed.population.dropped_count++; + digest = fnv1a64(ordered[i].second + "\n", digest); + if (ordered.size() - policy.max_labels_per_seed <= 1000) + seed.population.dropped.push_back(ordered[i].second); + } + } + seed.population.included = seed.labels.size(); + if (seed.population.dropped_count) + seed.population.dropped_digest = hex64(digest); + std::sort(seed.labels.begin(), seed.labels.end()); + seed.seed_id = make_seed_id(seed.sequence, canonical_orientation, seed.labels); + selection.seeds.push_back(std::move(seed)); + } + + for (size_t i = 0; i < selection.seeds.size(); ++i) { + for (size_t j = 0; j < selection.seeds.size(); ++j) { + const auto &a = selection.seeds[i].kmers, &b = selection.seeds[j].kmers; + if (i != j && a.begin < b.end && b.begin < a.end) + selection.seeds[i].overlaps_with.push_back(j); + } + } + return selection; +} + +} // namespace traversal +} // namespace graph +} // namespace mtg diff --git a/metagraph/src/graph/traversal/resolve.hpp b/metagraph/src/graph/traversal/resolve.hpp new file mode 100644 index 0000000000..9188997dc6 --- /dev/null +++ b/metagraph/src/graph/traversal/resolve.hpp @@ -0,0 +1,160 @@ +#ifndef __TRAVERSAL_RESOLVE_HPP__ +#define __TRAVERSAL_RESOLVE_HPP__ + +#include +#include +#include +#include +#include + +#include "traversal_types.hpp" +#include "label_oracle.hpp" + + +namespace mtg { +namespace graph { +namespace traversal { + +// half-open interval of k-mer starts [begin, end) +struct KmerInterval { + uint64_t begin = 0; + uint64_t end = 0; + + uint64_t size() const { return end - begin; } + bool operator==(const KmerInterval &o) const { return begin == o.begin && end == o.end; } + bool operator<(const KmerInterval &o) const { + return std::tie(begin, end) < std::tie(o.begin, o.end); + } +}; + +// The rows a discovery decodes at once +constexpr size_t kResolveBatchRows = 4096; + +struct ResolveOptions { + // exactly one of |labels| (explicit) or |discover| must be set + std::vector labels; + bool discover = false; + size_t discover_max_labels = 1000; + LabelKind discover_kind = LabelKind::COLUMN; + // the shortest support run (in k-mers) that becomes a seed candidate + uint64_t min_block_kmers = 1; +}; + +struct LabelProfile { + LabelRef label; + uint64_t kmers_supported = 0; + std::vector runs; // maximal support runs +}; + +struct SeedCandidate { + KmerInterval kmers; + std::vector labels; // indices into SupportProfile::labels, sorted +}; + +// A discovery that found more labels than |discover_max_labels| keeps those supporting the +// most k-mers (ties: smaller column, then sequence id) +struct LabelTruncation { + size_t kept = 0; + size_t total = 0; + uint64_t min_kept_kmers = 0; + uint64_t max_dropped_kmers = 0; + size_t dropped_full_length = 0; // dropped labels supporting every in-graph k-mer +}; + +struct SupportProfile { + size_t k = 0; + uint64_t num_kmers = 0; + Regime regime = Regime::BASIC; + std::vector graph_runs; // runs of k-mers present in the graph + std::vector labels; + std::optional labels_truncated; + std::vector candidates; // ordered: longer, more labels, smaller begin + + // base interval of a k-mer interval + std::pair bp_interval(const KmerInterval &iv) const { + return { iv.begin, iv.end + k - 1 }; + } +}; + +// Resolve where |query| is supported by which labels. No traversal. +// Throws std::invalid_argument for unknown labels or invalid options. +SupportProfile resolve_support(const LabelOracle &oracle, + std::string_view query, + const ResolveOptions &options); + +struct SelectionPolicy { + enum Policy { LONGEST_FIRST, MAX_SUPPORT }; + + Policy policy = MAX_SUPPORT; + size_t max_seeds = 10; + uint64_t min_block_bp = 0; // 0: at least k + size_t max_labels_per_seed = 1000; + bool merge_overlapping = true; +}; + +// The labels supporting a selected seed beyond |max_labels_per_seed|: the kept ones are the +// first in the order of their FNV-1a hash, so a cut is reproducible and unbiased by name +struct LabelPopulation { + size_t supporting_total = 0; // labels supporting the whole interval + size_t included = 0; + std::vector dropped; // kept only if <= 1000 + size_t dropped_count = 0; + std::string dropped_digest; +}; + +struct FrozenSeed { + std::string seed_id; + std::string sequence; + KmerInterval kmers; + std::vector labels; // bytewise sorted + LabelPopulation population; + std::vector overlaps_with; // indices of other seeds sharing k-mers +}; + +struct SeedSelection { + SelectionPolicy policy; + size_t num_candidates = 0; + size_t num_eligible = 0; + std::vector seeds; +}; + +SeedSelection select_seeds(const SupportProfile &profile, + std::string_view query, + const SelectionPolicy &policy, + bool canonical_orientation); + +// seed_id = hex FNV-1a-64 over the case-mapped (and in canonical graphs canonical) +// sequence and the length-prefixed sorted labels +std::string make_seed_id(std::string_view sequence, + bool canonical_orientation, + std::vector labels); + +uint64_t fnv1a64(std::string_view data, uint64_t hash = 0xcbf29ce484222325ULL); +// lower-case 16-digit hex of a 64-bit hash (seed ids, dropped-label digests) +std::string hex64(uint64_t x); + +// Run-length encoding of a presence mask ("xo..."), compatible with the +// query output of `with_signature`. +std::string encode_runs(const std::vector &runs, uint64_t num_kmers); + +// The maximal runs of the positions i in [0, n) where |present|(i) holds, ascending +template +std::vector runs_of(uint64_t n, const Present &present) { + std::vector runs; + for (uint64_t i = 0; i < n; ++i) { + if (!present(i)) + continue; + if (runs.empty() || runs.back().end != i) { + runs.push_back({ i, i + 1 }); + } else { + runs.back().end = i + 1; + } + } + return runs; +} + +} // namespace traversal +} // namespace graph +} // namespace mtg + +#endif // __TRAVERSAL_RESOLVE_HPP__ diff --git a/metagraph/src/graph/traversal/traversal_types.cpp b/metagraph/src/graph/traversal/traversal_types.cpp new file mode 100644 index 0000000000..089187b63d --- /dev/null +++ b/metagraph/src/graph/traversal/traversal_types.cpp @@ -0,0 +1,62 @@ +#include "traversal_types.hpp" + + +namespace mtg { +namespace graph { +namespace traversal { + +const char* to_string(Regime regime) { + switch (regime) { + case Regime::BASIC: return "basic"; + case Regime::PRIMARY: return "primary"; + case Regime::CANONICAL: return "canonical"; + } + return "unknown"; +} + +const char* to_string(LabelKind kind) { + switch (kind) { + case LabelKind::COLUMN: return "column"; + case LabelKind::HEADER: return "header"; + } + return "unknown"; +} + +const char* to_string(Arm arm) { + return arm == Arm::LEFT ? "left" : "right"; +} + +const char* to_string(EndReason reason) { + switch (reason) { + case EndReason::DEAD_END: return "dead_end"; + case EndReason::LABEL_LOST: return "label_lost"; + case EndReason::BRANCH: return "branch"; + case EndReason::EDGE_REUSE: return "edge_reuse"; + case EndReason::EDGE_REUSE_RC: return "edge_reuse_rc"; + case EndReason::REACHED_SEED: return "rejoined_seed"; + case EndReason::MAX_EXTENSION: return "max_extension_bp"; + case EndReason::MAX_STEPS: return "max_steps"; + case EndReason::MAX_LIVE_PATHS: return "max_live_paths"; + case EndReason::MAX_PATHS: return "max_paths"; + case EndReason::MAX_OUTPUT: return "max_output_bp"; + case EndReason::TIME_BUDGET: return "time_budget"; + } + return "unknown"; +} + +bool is_resource_stop(EndReason reason) { + switch (reason) { + case EndReason::MAX_STEPS: + case EndReason::MAX_LIVE_PATHS: + case EndReason::MAX_PATHS: + case EndReason::MAX_OUTPUT: + case EndReason::TIME_BUDGET: + return true; + default: + return false; + } +} + +} // namespace traversal +} // namespace graph +} // namespace mtg diff --git a/metagraph/src/graph/traversal/traversal_types.hpp b/metagraph/src/graph/traversal/traversal_types.hpp new file mode 100644 index 0000000000..505b0bf654 --- /dev/null +++ b/metagraph/src/graph/traversal/traversal_types.hpp @@ -0,0 +1,74 @@ +#ifndef __TRAVERSAL_TYPES_HPP__ +#define __TRAVERSAL_TYPES_HPP__ + +#include +#include + +#include "graph/representation/base/sequence_graph.hpp" +#include "annotation/binary_matrix/base/binary_matrix.hpp" + + +namespace mtg { +namespace graph { +namespace traversal { + +using node_index = SequenceGraph::node_index; +using Column = annot::matrix::BinaryMatrix::Column; +using Row = annot::matrix::BinaryMatrix::Row; +// index of a label in a per-seed label dictionary +using LabelId = uint32_t; + +static constexpr node_index npos = SequenceGraph::npos; + +// How oriented traversal nodes relate to annotation rows. +enum class Regime { + BASIC, // annotation key == node + PRIMARY, // CanonicalDBG wrapper over a primary graph: key = get_base_node(node) + CANONICAL, // native graph storing both strands: key = the canonical k-mer's node +}; +const char* to_string(Regime regime); + +// What a label refers to. +enum class LabelKind { + COLUMN, // an annotation column (e.g. a sample, file or taxid) + HEADER, // one sequence (FASTA header) inside a column, via CoordToHeader +}; +const char* to_string(LabelKind kind); + +struct LabelRef { + LabelKind kind = LabelKind::COLUMN; + Column column = 0; + uint64_t seq_id = 0; // only meaningful for HEADER + std::string name; // column name or FASTA header +}; + +enum class Arm : uint8_t { LEFT = 0, RIGHT = 1 }; +const char* to_string(Arm arm); + +// Why a label left a path, or why a path ended. +enum class EndReason : uint8_t { + // semantic, per label + DEAD_END = 0, // no structural successor in this direction + LABEL_LOST, // successors exist but none carries the label + BRANCH, // the label is on several successors and its branch limit is exhausted + EDGE_REUSE, // the label's continuation uses an edge already used on this path + EDGE_REUSE_RC, // ... used in the opposite orientation (hairpin / inverted repeat) + REACHED_SEED, // ... the edge leads back into the seed (e.g. a closed circle) + MAX_EXTENSION, // the requested extension length was reached + // caps: exploration of the requested extension is incomplete + MAX_STEPS, + MAX_LIVE_PATHS, + MAX_PATHS, + MAX_OUTPUT, + TIME_BUDGET, +}; +static constexpr size_t kNumEndReasons = static_cast(EndReason::TIME_BUDGET) + 1; +const char* to_string(EndReason reason); +// true if the reason means the requested extension was not explored completely +bool is_resource_stop(EndReason reason); + +} // namespace traversal +} // namespace graph +} // namespace mtg + +#endif // __TRAVERSAL_TYPES_HPP__ diff --git a/metagraph/tests/graph/traversal/test_label_oracle.cpp b/metagraph/tests/graph/traversal/test_label_oracle.cpp new file mode 100644 index 0000000000..a3c75d8054 --- /dev/null +++ b/metagraph/tests/graph/traversal/test_label_oracle.cpp @@ -0,0 +1,377 @@ +#include "gtest/gtest.h" + +#include +#include + +#include "tests/test_helpers.hpp" +#include "tests/graph/all/test_dbg_helpers.hpp" +#include "tests/annotation/test_annotated_dbg_helpers.hpp" + +#include "graph/traversal/label_oracle.hpp" +#include "graph/annotated_dbg.hpp" +#include "graph/representation/canonical_dbg.hpp" +#include "graph/representation/hash/dbg_sshash.hpp" +#include "annotation/coord_to_header.hpp" +#include "annotation/representation/column_compressed/annotate_column_compressed.hpp" +#include "annotation/representation/annotation_matrix/static_annotators_def.hpp" +#include "common/seq_tools/reverse_complement.hpp" + + +namespace { + +using namespace mtg; +using namespace mtg::graph; +using namespace mtg::graph::traversal; +using namespace mtg::test; + +const size_t kTestK = 11; + +// Fixed random sequences: long enough for every k-mer to be distinct with high +// probability, short enough for RowDiff fixtures to build quickly. +std::vector make_sequences(size_t n, size_t len, uint32_t seed) { + std::mt19937 gen(seed); + std::vector seqs; + for (size_t i = 0; i < n; ++i) { + std::string s(len, 'A'); + for (char &c : s) { + c = "ACGT"[gen() % 4]; + } + seqs.push_back(s); + } + return seqs; +} + +#if ! _PROTEIN_GRAPH +std::string rc(std::string s) { ::reverse_complement(s); return s; } +#endif + +std::vector all_modes() { + return { + DeBruijnGraph::BASIC, +#if ! _PROTEIN_GRAPH + DeBruijnGraph::CANONICAL, + DeBruijnGraph::PRIMARY, +#endif + }; +} + +template +class LabelOracleTest : public ::testing::Test {}; + +typedef ::testing::Types< + std::pair>, + std::pair>, + std::pair>, + std::pair>, + std::pair>, + std::pair, + std::pair, + std::pair, + std::pair, + std::pair +> OracleTypes; +TYPED_TEST_SUITE(LabelOracleTest, OracleTypes); + +Regime expected_regime(const DeBruijnGraph &graph, DeBruijnGraph::Mode mode) { + if (dynamic_cast(&graph)) + return Regime::PRIMARY; + // DBGSSHash builds PRIMARY as native CANONICAL and is not wrapped + if (graph.get_mode() == DeBruijnGraph::CANONICAL) + return Regime::CANONICAL; + EXPECT_EQ(DeBruijnGraph::BASIC, mode); + return Regime::BASIC; +} + +// regime detection and oriented-node -> annotation-key mapping +TYPED_TEST(LabelOracleTest, RegimeAndKeyMapping) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + auto seqs = make_sequences(3, 80, 1); + std::vector labels { "A", "B", "C" }; + + for (DeBruijnGraph::Mode mode : all_modes()) { + auto anno = build_anno_graph(kTestK, seqs, labels, mode); + LabelOracle oracle(*anno); + EXPECT_EQ(expected_regime(anno->get_graph(), mode), oracle.regime()); + EXPECT_EQ(kTestK, oracle.get_k()); + EXPECT_EQ(3u, oracle.num_columns()); + + const DeBruijnGraph &graph = oracle.graph(); + for (size_t s = 0; s < seqs.size(); ++s) { + const std::string &seq = seqs[s]; + auto keys = oracle.keys_of_sequence(seq); + auto nodes = map_to_nodes_sequentially(graph, seq); + ASSERT_EQ(seq.size() - kTestK + 1, keys.size()); + ASSERT_EQ(keys.size(), nodes.size()); + + LabelQuery query(oracle, { oracle.resolve_label(labels[s]) }); + EXPECT_STREQ("rows", query.access_path()); + auto hits = query.fetch(keys); + for (size_t i = 0; i < keys.size(); ++i) { + ASSERT_NE(npos, nodes[i]); + ASSERT_NE(npos, keys[i]); + std::string_view kmer(seq.data() + i, kTestK); + EXPECT_EQ(keys[i], oracle.key_of(nodes[i], kmer)); + // the label written on this k-mer is visible through the key + EXPECT_EQ((LabelQuery::NodeHits{ 0 }), hits[i]) << "k-mer " << kmer; + } + +#if ! _PROTEIN_GRAPH + // reverse-complement orientation + std::string seq_rc = rc(seq); + auto nodes_rc = map_to_nodes_sequentially(graph, seq_rc); + auto keys_rc = oracle.keys_of_sequence(seq_rc); + std::reverse(nodes_rc.begin(), nodes_rc.end()); + std::reverse(keys_rc.begin(), keys_rc.end()); + auto hits_rc = query.fetch(keys_rc); + if (oracle.regime() == Regime::BASIC) { + // forward strand only: an RC k-mer is in the graph only if it occurs on the + // forward strand of some sequence, and then it carries that sequence's label + for (size_t i = 0; i < keys_rc.size(); ++i) { + std::string rc_kmer = rc(seq.substr(i, kTestK)); + bool on_forward_strand = false; + for (const auto &other : seqs) { + on_forward_strand |= other.find(rc_kmer) != std::string::npos; + } + EXPECT_EQ(on_forward_strand, keys_rc[i] != npos) << "rc k-mer " << rc_kmer; + bool in_this_seq = seq.find(rc_kmer) != std::string::npos; + EXPECT_EQ(in_this_seq, !hits_rc[i].empty()) << "rc k-mer " << rc_kmer; + } + } else { + // both orientations map to the same annotation key and see the label + for (size_t i = 0; i < keys.size(); ++i) { + ASSERT_NE(npos, nodes_rc[i]); + EXPECT_EQ(keys[i], keys_rc[i]); + std::string_view kmer_rc(seq_rc.data() + seq_rc.size() - kTestK - i, kTestK); + EXPECT_EQ(keys[i], oracle.key_of(nodes_rc[i], kmer_rc)); + EXPECT_EQ((LabelQuery::NodeHits{ 0 }), hits_rc[i]); + } + } +#endif + } + } +} + +// the hits agree with AnnotatedDBG::get_labels for every label subset +TYPED_TEST(LabelOracleTest, HitsAgreeWithGetLabels) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + // overlapping sequences so that k-mers carry several labels + auto base = make_sequences(1, 120, 7)[0]; + std::vector seqs { base, base.substr(20, 70), base.substr(50) + "ACGTACGTACGTAAACCC", + make_sequences(1, 60, 9)[0] }; + std::vector labels { "L0", "L1", "L2", "L3" }; + + for (DeBruijnGraph::Mode mode : all_modes()) { + auto anno = build_anno_graph(kTestK, seqs, labels, mode); + LabelOracle oracle(*anno); + + std::vector> label_sets { + { "L0" }, { "L1", "L3" }, { "L0", "L1", "L2", "L3" } + }; + for (const auto &names : label_sets) { + std::vector refs; + for (const auto &n : names) { + refs.push_back(oracle.resolve_label(n)); + } + LabelQuery query(oracle, refs); + for (const auto &seq : seqs) { + auto keys = oracle.keys_of_sequence(seq); + auto hits = query.fetch(keys); + ASSERT_EQ(keys.size(), hits.size()); + for (size_t i = 0; i < keys.size(); ++i) { + std::vector expected; + for (const auto &l : anno->get_labels(seq.substr(i, kTestK), 1.0)) { + if (std::find(names.begin(), names.end(), l) != names.end()) + expected.push_back(l); + } + std::vector got; + for (LabelId h : hits[i]) { + got.push_back(names[h]); + } + std::sort(expected.begin(), expected.end()); + std::sort(got.begin(), got.end()); + EXPECT_EQ(expected, got) << "k-mer " << seq.substr(i, kTestK); + } + } + } + EXPECT_GT(oracle.counters().rows_requested, 0u); + EXPECT_GT(oracle.counters().cache_hits, 0u); + EXPECT_THROW(oracle.resolve_label("no_such_label"), std::invalid_argument); + EXPECT_THROW(LabelQuery(oracle, {}), std::invalid_argument); + EXPECT_THROW(LabelQuery(oracle, { oracle.resolve_label("L0"), oracle.resolve_label("L0") }), + std::invalid_argument); + } +} + +// Unmasked DBGSuccinct (CLI-built graphs keep dummy k-mers): keys never point at dummies +TEST(LabelOracle, UnmaskedSuccinctKeys) { + auto seqs = make_sequences(2, 60, 3); + for (DeBruijnGraph::Mode mode : { DeBruijnGraph::BASIC, DeBruijnGraph::CANONICAL }) { + auto graph = std::make_shared(kTestK, mode); + for (const auto &s : seqs) { + graph->add_sequence(s); + } + ASSERT_EQ(nullptr, graph->get_mask()); + auto anno = std::make_unique( + graph, std::make_unique>(graph->max_index())); + anno->annotate_sequence(seqs[0], { "A" }); + anno->annotate_sequence(seqs[1], { "B" }); + + LabelOracle oracle(*anno); + EXPECT_EQ(mode == DeBruijnGraph::BASIC ? Regime::BASIC : Regime::CANONICAL, oracle.regime()); + LabelQuery query(oracle, { oracle.resolve_label("A"), oracle.resolve_label("B") }); + auto keys = oracle.keys_of_sequence(seqs[0]); + auto hits = query.fetch(keys); + for (size_t i = 0; i < keys.size(); ++i) { + ASSERT_NE(npos, keys[i]); + EXPECT_LE(keys[i], graph->max_index()); + ASSERT_GE(hits[i].size(), 1u); + EXPECT_EQ(0u, hits[i][0]); + } + } +} + +// A batch that mixes cached and uncached keys must survive eviction: the cache is cleared +// when it would overflow, so every key of the current call has to be refetched +TEST(LabelOracle, EvictionKeepsTheCurrentBatchAnswerable) { + auto seqs = make_sequences(1, 60, 21); + auto anno = build_anno_graph>( + kTestK, seqs, { "A" }, DeBruijnGraph::BASIC); + LabelOracle oracle(*anno); + auto keys = oracle.keys_of_sequence(seqs[0]); + ASSERT_GE(keys.size(), 3u); + + for (size_t capacity : { size_t(1), size_t(2), size_t(3) }) { + LabelQuery query(oracle, { oracle.resolve_label("A") }, capacity); + auto first = query.fetch({ keys[0], keys[1] }); + ASSERT_EQ(2u, first.size()); + // keys[0] is a hit, keys[2] a miss: the mixed batch forces an eviction + auto second = query.fetch({ keys[0], keys[2] }); + ASSERT_EQ(2u, second.size()); + EXPECT_EQ(first[0], second[0]) << "capacity " << capacity; + for (const auto &hits : second) { + EXPECT_EQ((LabelQuery::NodeHits{ 0 }), hits); + } + // and a batch larger than the whole cache is still answered + auto all = query.fetch(keys); + ASSERT_EQ(keys.size(), all.size()); + for (const auto &hits : all) { + EXPECT_EQ((LabelQuery::NodeHits{ 0 }), hits); + } + } +} + + +// Header labels via CoordToHeader (the refseq33m shape) +template +class LabelOracleCoordTest : public ::testing::Test {}; + +typedef ::testing::Types< + std::pair>, + std::pair>, + std::pair, + std::pair +> CoordOracleTypes; +TYPED_TEST_SUITE(LabelOracleCoordTest, CoordOracleTypes); + +TYPED_TEST(LabelOracleCoordTest, HeaderLabels) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + // one column "F" holding two sequences, like a FASTA file with two records, + // plus a second column "G" with one sequence that shares a 30 bp block with s1 + auto rnd = make_sequences(3, 70, 11); + std::string s1 = rnd[0], s2 = rnd[1], s3 = rnd[2].substr(0, 20) + s1.substr(10, 30) + rnd[2].substr(50); + std::vector seqs { s1, s2, s3 }; + std::vector labels { "F", "F", "G" }; + uint64_t n1 = s1.size() - kTestK + 1, n2 = s2.size() - kTestK + 1, n3 = s3.size() - kTestK + 1; + + auto anno = build_anno_graph(kTestK, seqs, labels, DeBruijnGraph::BASIC, + true, { 0, n1, 0 }); + ASSERT_TRUE(dynamic_cast(&anno->get_annotator().get_matrix())); + + // CoordToHeader: column order is the label encoder order + Column col_f = anno->get_annotator().get_label_encoder().encode("F"); + Column col_g = anno->get_annotator().get_label_encoder().encode("G"); + std::vector> headers(2); + std::vector> num_kmers(2); + headers[col_f] = { "acc1", "acc2" }; num_kmers[col_f] = { n1, n2 }; + headers[col_g] = { "acc3" }; num_kmers[col_g] = { n3 }; + annot::CoordToHeader cth(std::move(headers), std::move(num_kmers)); + + LabelOracle oracle(*anno, &cth); + ASSERT_TRUE(oracle.has_coordinates()); + ASSERT_EQ(&cth, oracle.coord_to_header()); + + auto acc1 = oracle.resolve_label("acc1"); + auto acc2 = oracle.resolve_label("acc2"); + auto acc3 = oracle.resolve_label("acc3"); + auto col_label_f = oracle.resolve_label("F"); + EXPECT_EQ(LabelKind::HEADER, acc1.kind); + EXPECT_EQ(col_f, acc1.column); + EXPECT_EQ(0u, acc1.seq_id); + EXPECT_EQ(1u, acc2.seq_id); + EXPECT_EQ(LabelKind::COLUMN, col_label_f.kind); + EXPECT_EQ(std::make_pair(n1, n1 + n2 - 1), oracle.sequence_coords(col_f, 1)); + EXPECT_EQ(std::make_pair(uint64_t(1), uint64_t(3)), oracle.map_coord(col_f, n1 + 3)); + + LabelQuery query(oracle, { acc1, acc2, acc3, col_label_f }); + EXPECT_STREQ("tuples", query.access_path()); + + // s2's k-mers: header acc2 and column F, no acc1 + auto hits2 = query.fetch(oracle.keys_of_sequence(s2)); + for (size_t i = 0; i < hits2.size(); ++i) { + EXPECT_EQ((LabelQuery::NodeHits{ 1, 3 }), hits2[i]) << "s2 k-mer " << i; + } + + // s1's k-mers: acc1 everywhere; acc3 on the k-mers shared with s3 + auto hits1 = query.fetch(oracle.keys_of_sequence(s1)); + for (size_t i = 0; i < hits1.size(); ++i) { + bool shared = s3.find(s1.substr(i, kTestK)) != std::string::npos; + const LabelQuery::NodeHits expected = shared ? LabelQuery::NodeHits{ 0, 2, 3 } + : LabelQuery::NodeHits{ 0, 3 }; + EXPECT_EQ(expected, hits1[i]) << "s1 k-mer " << i; + } + EXPECT_GT(oracle.counters().tuple_rows_fetched, 0u); + + // header labels need a CoordToHeader + LabelOracle no_headers(*anno); + EXPECT_THROW(no_headers.resolve_label("acc1"), std::invalid_argument); + EXPECT_THROW(LabelQuery(no_headers, { acc1 }), std::invalid_argument); +} + +// The reverse header index answers each header's first record, also after a load() +TEST(CoordToHeaderIndex, FindHeader) { + auto anno = build_anno_graph>( + 5, { "ACGTAGATCGAA", "GTCTAACGTTGA", "TATGCCAGATCA", "AGCGTCTTGCGA" }, + { "F", "F", "F", "F" }, DeBruijnGraph::BASIC, true, { 0, 8, 16, 24 }); + for (const std::vector &order : { std::vector{ "A", "X", "Y", "Z" }, + std::vector{ "A", "Z", "Y", "X" } }) { + annot::CoordToHeader cth({ order }, { { 8, 8, 8, 8 } }); + LabelOracle oracle(*anno, &cth); + for (size_t s = 0; s < order.size(); ++s) { + auto ref = oracle.find_header(order[s]); + ASSERT_TRUE(ref) << order[s]; + EXPECT_EQ(s, ref->seq_id) << order[s]; + EXPECT_EQ(order[s], oracle.header_name(ref->column, ref->seq_id)); + EXPECT_EQ(ref->seq_id, oracle.resolve_label(order[s]).seq_id); + } + EXPECT_FALSE(oracle.find_header("W")); + } + + // a header repeated in a second column resolves to its first column + annot::CoordToHeader two({ { "A", "B" }, { "B", "C" } }, { { 1, 1 }, { 1, 1 } }); + EXPECT_EQ(std::make_pair(Column(0), size_t(1)), *two.find_header("B")); + EXPECT_EQ(std::make_pair(Column(1), size_t(1)), *two.find_header("C")); + + // load() replaces the headers and the index with them + const std::string base = test_dump_dir() + "/coord_to_header_index"; + two.serialize(base); + annot::CoordToHeader loaded({ { "Q" } }, { { 1 } }); + ASSERT_TRUE(loaded.find_header("Q")); + ASSERT_TRUE(loaded.load(base)); + EXPECT_FALSE(loaded.find_header("Q")); + EXPECT_EQ(std::make_pair(Column(1), size_t(1)), *loaded.find_header("C")); +} + +} // namespace diff --git a/metagraph/tests/graph/traversal/test_resolve.cpp b/metagraph/tests/graph/traversal/test_resolve.cpp new file mode 100644 index 0000000000..82a235ec53 --- /dev/null +++ b/metagraph/tests/graph/traversal/test_resolve.cpp @@ -0,0 +1,403 @@ +#include "gtest/gtest.h" + +#include +#include +#include +#include + +#include "tests/graph/all/test_dbg_helpers.hpp" +#include "tests/annotation/test_annotated_dbg_helpers.hpp" + +#include "graph/traversal/resolve.hpp" +#include "graph/annotated_dbg.hpp" +#include "annotation/coord_to_header.hpp" +#include "annotation/representation/column_compressed/annotate_column_compressed.hpp" +#include "annotation/representation/annotation_matrix/static_annotators_def.hpp" + + +namespace { + +using namespace mtg; +using namespace mtg::graph; +using namespace mtg::graph::traversal; +using namespace mtg::test; + +const size_t kK = 11; + +std::string random_seq(size_t len, uint32_t seed) { + std::mt19937 gen(seed); + std::string s(len, 'A'); + for (char &c : s) { + c = "ACGT"[gen() % 4]; + } + return s; +} + +std::vector iv(std::initializer_list> l) { + std::vector out; + for (auto [a, b] : l) { + out.push_back({ a, b }); + } + return out; +} + +template +class ResolveTest : public ::testing::Test {}; + +typedef ::testing::Types< + std::pair>, + std::pair>, + std::pair>, + std::pair, + std::pair +> ResolveTypes; +TYPED_TEST_SUITE(ResolveTest, ResolveTypes); + +// Query = 200 bp with a 15 bp foreign insert at [80,95); label A is written on +// query[0,120) without the insert, B on query[100,200). +TYPED_TEST(ResolveTest, ThreeStatesAndRuns) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + std::string q = random_seq(200, 42); + std::string a_graph = q.substr(0, 80) + q.substr(95, 25); + std::string b_seq = q.substr(100); + + for (DeBruijnGraph::Mode mode : { DeBruijnGraph::BASIC +#if ! _PROTEIN_GRAPH + , DeBruijnGraph::CANONICAL, DeBruijnGraph::PRIMARY +#endif + }) { + auto anno = build_anno_graph(kK, { a_graph, b_seq }, { "A", "B" }, mode); + LabelOracle oracle(*anno); + const bool canonical = mode != DeBruijnGraph::BASIC; + + ResolveOptions opts; + opts.labels = { "A", "B" }; + auto profile = resolve_support(oracle, q, opts); + EXPECT_EQ(190u, profile.num_kmers); + EXPECT_EQ(kK, profile.k); + + // k-mers overlapping the insert [80,95) are absent from the graph: starts 70..94 + ASSERT_EQ(2u, profile.graph_runs.size()) << "mode " << mode; + EXPECT_EQ((KmerInterval{ 0, 70 }), profile.graph_runs[0]); + EXPECT_EQ((KmerInterval{ 95, 190 }), profile.graph_runs[1]); + + ASSERT_EQ(2u, profile.labels.size()); + const auto &A = profile.labels[0], &B = profile.labels[1]; + EXPECT_EQ("A", A.label.name); + // A: present on [0,70) and on [95,110) (k-mers fully inside q[95,120)) + EXPECT_EQ(iv({ {0, 70}, {95, 110} }), A.runs); + EXPECT_EQ(85u, A.kmers_supported); + // B: k-mers fully inside q[100,200): starts 100..189 + EXPECT_EQ(iv({ {100, 190} }), B.runs); + EXPECT_EQ(90u, B.kmers_supported); + + // candidates: identical runs grouped, ordered longer first + ASSERT_EQ(3u, profile.candidates.size()); + EXPECT_EQ((KmerInterval{ 100, 190 }), profile.candidates[0].kmers); + EXPECT_EQ((std::vector{ 1 }), profile.candidates[0].labels); + EXPECT_EQ((KmerInterval{ 0, 70 }), profile.candidates[1].kmers); + EXPECT_EQ((KmerInterval{ 95, 110 }), profile.candidates[2].kmers); + + EXPECT_EQ("x70o25x15o80", encode_runs(A.runs, profile.num_kmers)); + + // longest_first + SelectionPolicy pol; + pol.policy = SelectionPolicy::LONGEST_FIRST; + pol.max_seeds = 2; + pol.merge_overlapping = false; + auto sel = select_seeds(profile, q, pol, canonical); + EXPECT_EQ(3u, sel.num_candidates); + EXPECT_EQ(3u, sel.num_eligible); + ASSERT_EQ(2u, sel.seeds.size()); + EXPECT_EQ((KmerInterval{ 100, 190 }), sel.seeds[0].kmers); + EXPECT_EQ(q.substr(100), sel.seeds[0].sequence); + EXPECT_EQ((std::vector{ "B" }), sel.seeds[0].labels); + EXPECT_EQ((KmerInterval{ 0, 70 }), sel.seeds[1].kmers); + EXPECT_EQ(q.substr(0, 80), sel.seeds[1].sequence); + EXPECT_TRUE(sel.seeds[0].overlaps_with.empty()); + EXPECT_EQ(16u, sel.seeds[0].seed_id.size()); + EXPECT_NE(sel.seeds[0].seed_id, sel.seeds[1].seed_id); + + // max_support with overlap: A and B both cover [100,110) + pol = SelectionPolicy(); + pol.max_seeds = 1; + sel = select_seeds(profile, q, pol, canonical); + ASSERT_EQ(1u, sel.seeds.size()); + EXPECT_EQ((KmerInterval{ 100, 110 }), sel.seeds[0].kmers); + EXPECT_EQ((std::vector{ "A", "B" }), sel.seeds[0].labels); + EXPECT_EQ(2u, sel.seeds[0].population.supporting_total); + + // max_labels_per_seed drops deterministically and reports the drop + pol.max_labels_per_seed = 1; + sel = select_seeds(profile, q, pol, canonical); + ASSERT_EQ(1u, sel.seeds.size()); + EXPECT_EQ(1u, sel.seeds[0].labels.size()); + EXPECT_EQ(1u, sel.seeds[0].population.dropped_count); + EXPECT_EQ(1u, sel.seeds[0].population.dropped.size()); + EXPECT_FALSE(sel.seeds[0].population.dropped_digest.empty()); + EXPECT_EQ(sel.seeds[0].labels, select_seeds(profile, q, pol, canonical).seeds[0].labels); + + // agreement with the signature API for the labels present + auto sigs = anno->get_top_label_signatures(q, 10, 0.0, 0.0); + for (const auto &[label, count, mask] : sigs) { + const auto &lp = label == "A" ? A : B; + EXPECT_EQ(count, lp.kmers_supported); + std::vector expected(mask.size()); + for (size_t i = 0; i < mask.size(); ++i) { + expected[i] = mask[i]; + } + std::vector got(profile.num_kmers, false); + for (const auto &run : lp.runs) { + for (uint64_t i = run.begin; i < run.end; ++i) { + got[i] = true; + } + } + EXPECT_EQ(expected, got) << label; + } + } +} + +// discovery with truncation keeps the top labels by k-mers, ties by column id +TYPED_TEST(ResolveTest, DiscoverTruncation) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + std::string q = random_seq(100, 5); + // X, Y, Z carry the whole query; W carries half + auto anno = build_anno_graph(kK, { q, q, q, q.substr(0, 50) }, + { "X", "Y", "Z", "W" }, DeBruijnGraph::BASIC); + LabelOracle oracle(*anno); + ResolveOptions opts; + opts.discover = true; + opts.discover_max_labels = 2; + auto profile = resolve_support(oracle, q, opts); + ASSERT_TRUE(profile.labels_truncated.has_value()); + EXPECT_EQ(2u, profile.labels_truncated->kept); + EXPECT_EQ(4u, profile.labels_truncated->total); + EXPECT_EQ(90u, profile.labels_truncated->min_kept_kmers); + EXPECT_EQ(90u, profile.labels_truncated->max_dropped_kmers); + EXPECT_EQ(1u, profile.labels_truncated->dropped_full_length); // Z dropped, W is not full-length + ASSERT_EQ(2u, profile.labels.size()); + Column cx = anno->get_annotator().get_label_encoder().encode("X"); + Column cy = anno->get_annotator().get_label_encoder().encode("Y"); + EXPECT_EQ(std::min(cx, cy), profile.labels[0].label.column); + EXPECT_EQ(std::max(cx, cy), profile.labels[1].label.column); + + opts.discover_max_labels = 10; + profile = resolve_support(oracle, q, opts); + EXPECT_FALSE(profile.labels_truncated.has_value()); + ASSERT_EQ(4u, profile.labels.size()); + EXPECT_EQ("W", profile.labels[3].label.name); + EXPECT_EQ(40u, profile.labels[3].kmers_supported); + + // validation + ResolveOptions bad; + EXPECT_THROW(resolve_support(oracle, q, bad), std::invalid_argument); + bad.labels = { "nope" }; + EXPECT_THROW(resolve_support(oracle, q, bad), std::invalid_argument); + bad.discover = true; + EXPECT_THROW(resolve_support(oracle, q, bad), std::invalid_argument); + ResolveOptions headers; + headers.discover = true; + headers.discover_kind = LabelKind::HEADER; + EXPECT_THROW(resolve_support(oracle, q, headers), std::invalid_argument); // no CoordToHeader + EXPECT_THROW(resolve_support(oracle, q.substr(0, kK - 1), opts), std::invalid_argument); +} + +// max_support at many labels: a synthetic profile whose runs end at jittered positions. The +// breadth-optimal block is the prefix every label covers +TEST(Resolve, MaxSupportScalesWithManyLabels) { + const size_t num_labels = 1000; + const uint64_t num_kmers = 5000; + SupportProfile profile; + profile.k = kK; + profile.num_kmers = num_kmers; + profile.graph_runs = iv({ {0, num_kmers} }); + std::mt19937 gen(7); + for (size_t i = 0; i < num_labels; ++i) { + LabelProfile lp; + lp.label.column = static_cast(i); + lp.label.name = "L" + std::to_string(i); + // every label covers [0, 1000); each ends somewhere in [1000, 5000) + uint64_t end = 1000 + gen() % (num_kmers - 1000); + lp.runs = iv({ {0, end} }); + lp.kmers_supported = end; + profile.labels.push_back(lp); + } + for (LabelId l = 0; l < profile.labels.size(); ++l) { + profile.candidates.push_back({ profile.labels[l].runs[0], { l } }); + } + + SelectionPolicy policy; + policy.max_seeds = 1; + policy.min_block_bp = kK; // min_kmers == 1 + const std::string query(num_kmers + kK - 1, 'A'); + auto selection = select_seeds(profile, query, policy, false); + ASSERT_EQ(1u, selection.seeds.size()); + uint64_t smallest_end = num_kmers; + for (const auto &lp : profile.labels) { + smallest_end = std::min(smallest_end, lp.runs[0].end); + } + EXPECT_EQ((KmerInterval{ 0, smallest_end }), selection.seeds[0].kmers); + EXPECT_EQ(num_labels, selection.seeds[0].labels.size()); + + // a length floor trades labels for length + policy.min_block_bp = 3000 + kK - 1; + auto longer = select_seeds(profile, query, policy, false); + ASSERT_EQ(1u, longer.seeds.size()); + EXPECT_GE(longer.seeds[0].kmers.size(), 3000u); + EXPECT_LT(longer.seeds[0].labels.size(), num_labels); + for (const auto &name : longer.seeds[0].labels) { + const auto &lp = *std::find_if(profile.labels.begin(), profile.labels.end(), + [&](const LabelProfile &x) { return x.label.name == name; }); + EXPECT_GE(lp.runs[0].end, longer.seeds[0].kmers.end); + } +} + +// Later max_support rounds see the runs CLIPPED by what earlier rounds took: with A on +// [0,10) and B on [0,5), the first round takes [0,5) (2 labels), the second [5,10) +TEST(Resolve, MaxSupportClipsRunsBetweenRounds) { + SupportProfile profile; + profile.k = kK; + profile.num_kmers = 10; + profile.graph_runs = iv({ {0, 10} }); + for (auto [name, end] : { std::pair{ "A", 10 }, + std::pair{ "B", 5 } }) { + LabelProfile lp; + lp.label.column = static_cast(profile.labels.size()); + lp.label.name = name; + lp.runs = iv({ {0, end} }); + lp.kmers_supported = end; + profile.labels.push_back(lp); + } + + SelectionPolicy policy; + policy.max_seeds = 2; + policy.min_block_bp = kK; // min_kmers == 1 + policy.merge_overlapping = false; // keep the rounds' intervals as chosen + auto sel = select_seeds(profile, std::string(10 + kK - 1, 'A'), policy, false); + ASSERT_EQ(2u, sel.seeds.size()); + EXPECT_EQ((KmerInterval{ 0, 5 }), sel.seeds[0].kmers); + EXPECT_EQ((std::vector{ "A", "B" }), sel.seeds[0].labels); + EXPECT_EQ((KmerInterval{ 5, 10 }), sel.seeds[1].kmers); + EXPECT_EQ((std::vector{ "A" }), sel.seeds[1].labels); + EXPECT_TRUE(sel.seeds[0].overlaps_with.empty()); + EXPECT_TRUE(sel.seeds[1].overlaps_with.empty()); +} + +// seed ids: canonical orientation and label order do not change the id +TEST(Resolve, SeedId) { + std::string s = "ACGTTGCAAGT"; + std::string rc = "ACTTGCAACGT"; + EXPECT_EQ(make_seed_id(s, true, { "b", "a" }), make_seed_id(rc, true, { "a", "b" })); + EXPECT_NE(make_seed_id(s, false, { "a" }), make_seed_id(rc, false, { "a" })); + // length-prefixed labels: {"a,b"} != {"a","b"} + EXPECT_NE(make_seed_id(s, false, { "a,b" }), make_seed_id(s, false, { "a", "b" })); + // the id of the labeled-traversal prototype for the same seed (no release name) + EXPECT_EQ(hex64(fnv1a64("1:a", fnv1a64("\t", fnv1a64(s, fnv1a64("\t"))))), + make_seed_id(s, false, { "a" })); + EXPECT_EQ("x3o2x5", encode_runs(iv({ {0, 3}, {5, 10} }), 10)); +} + + +template +class ResolveCoordTest : public ::testing::Test {}; +typedef ::testing::Types< + std::pair>, + std::pair +> ResolveCoordTypes; +TYPED_TEST_SUITE(ResolveCoordTest, ResolveCoordTypes); + +// header labels: acc1 = L + R + M + R (repeat R twice) and acc2 = R, both in column F +TYPED_TEST(ResolveCoordTest, HeaderLabelsAndDiscovery) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + std::string L = random_seq(40, 1), R = random_seq(40, 2), M = random_seq(40, 3); + std::string acc1 = L + R + M + R; + std::string acc2 = R; + uint64_t n1 = acc1.size() - kK + 1, n2 = acc2.size() - kK + 1; + auto anno = build_anno_graph(kK, { acc1, acc2 }, { "F", "F" }, + DeBruijnGraph::BASIC, true, { 0, n1 }); + annot::CoordToHeader cth({ { "acc1", "acc2" } }, { { n1, n2 } }); + LabelOracle oracle(*anno, &cth); + + std::string q = R + M + R; + ResolveOptions opts; + opts.labels = { "acc1", "acc2", "F" }; + auto profile = resolve_support(oracle, q, opts); + ASSERT_EQ(3u, profile.labels.size()); + EXPECT_TRUE(profile.labels[0].label.kind == LabelKind::HEADER); + EXPECT_EQ(iv({ {0, q.size() - kK + 1} }), profile.labels[0].runs); // acc1 + // acc2 = R: k-mers of the two R copies, i.e. [0,30) and [80,110) + EXPECT_EQ(iv({ {0, 30}, {80, 110} }), profile.labels[1].runs); + EXPECT_TRUE(profile.labels[2].label.kind == LabelKind::COLUMN); + EXPECT_EQ(iv({ {0, q.size() - kK + 1} }), profile.labels[2].runs); // column F + + // header discovery gives the same profiles + ResolveOptions disc; + disc.discover = true; + disc.discover_kind = LabelKind::HEADER; + auto discovered = resolve_support(oracle, q, disc); + ASSERT_EQ(2u, discovered.labels.size()); + EXPECT_EQ("acc1", discovered.labels[0].label.name); + EXPECT_EQ(profile.labels[0].runs, discovered.labels[0].runs); + EXPECT_EQ("acc2", discovered.labels[1].label.name); + EXPECT_EQ(60u, discovered.labels[1].kmers_supported); + EXPECT_EQ(profile.labels[1].runs, discovered.labels[1].runs); + + // an unknown header + opts.labels = { "acc3" }; + EXPECT_THROW(resolve_support(oracle, q, opts), std::invalid_argument); +} + +// The explicit profiles of a few hundred labels are the discovery's, and a label absent from +// the query has none +TYPED_TEST(ResolveCoordTest, ExplicitLabelsGiveTheDiscoverysProfiles) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + const std::string R = random_seq(40, 21), M = random_seq(40, 22); + const std::string q = random_seq(150, 23) + R + M + R + M + random_seq(150, 24); + std::mt19937 gen(25); + std::vector seqs, labels; + for (size_t i = 0; i < 300; ++i) { + // a stretch of the query, some with a second stretch after it, some with a base changed + const size_t a = gen() % (q.size() - 40); + std::string s = q.substr(a, 20 + gen() % std::min(200, q.size() - a - 20)); + if (i % 3 == 0) { + const size_t c = gen() % (q.size() - 40); + s += q.substr(c, 20 + gen() % 20); + } + if (i % 5 == 0) + s[gen() % s.size()] = "ACGT"[gen() % 4]; + seqs.push_back(s); + labels.push_back("L" + std::to_string(i)); + } + seqs.push_back(random_seq(60, 26)); + labels.push_back("absent"); + auto anno = build_anno_graph(kK, seqs, labels, DeBruijnGraph::BASIC, true); + LabelOracle oracle(*anno); + + ResolveOptions disc; + disc.discover = true; + const SupportProfile discovered = resolve_support(oracle, q, disc); + ResolveOptions opts; + opts.labels = labels; + const SupportProfile explicit_ = resolve_support(oracle, q, opts); + ASSERT_EQ(labels.size(), explicit_.labels.size()); + std::map by_name; + for (const LabelProfile &lp : explicit_.labels) { + by_name[lp.label.name] = &lp; + } + ASSERT_GT(discovered.labels.size(), 200u); + for (const LabelProfile &d : discovered.labels) { + const LabelProfile &e = *by_name.at(d.label.name); + EXPECT_EQ(d.runs, e.runs) << d.label.name; + EXPECT_EQ(d.kmers_supported, e.kmers_supported) << d.label.name; + } + const LabelProfile &absent = *by_name.at("absent"); + EXPECT_TRUE(absent.runs.empty()); + EXPECT_EQ(0u, absent.kmers_supported); + EXPECT_EQ(discovered.candidates.size(), explicit_.candidates.size()); +} + +} // namespace From 9daa8dae62ceff9bf1ef22e97bf3b1152ef524e1 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sun, 1 Nov 2026 23:12:26 +0100 Subject: [PATCH 02/52] Add the labeled traversal walker (v1) MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Second piece of the labeled-traversal v1 stack: extend a seed under a permitted label set, rebuilt from Gunnar's walker on gr/labeled-traversal (d5d26baa) for the v1 feature set only: - the seed's labels named explicitly (labels not carrying every seed k-mer are dropped and reported) or derived as the labels on every seed k-mer, column or header, capped with a dropped count and digest; - left, right or both arms, level-synchronous; - per path the labels carrying every k-mer walked; a label ends with dead_end, label_lost, branch (over its branch allowance), edge_reuse, edge_reuse_rc or rejoined_seed; hairpins are skipped; - a tree of segments per arm (paths are kept apart where they meet); - the extension length, the time budget and fixed safety caps (steps, live paths, paths, output bases, splits per path), with complete_to_bp stating where every admissible walk is present. Not carried over: label switching and loss budgets, trace support and record coordinates, annotate mode, memory and work budgets, read pacing, the lookahead, the attempt protocol, reconvergence merging, beams, quorums, continuations, growth profiles and branch-event refusals. Tests: the walker suite trimmed to the v1 features, and a check of every walker claim against an independent string-level model of the records (test_trie_reference.hpp: label_trie, from Gunnar's trie reference). Source paths on d5d26baa: metagraph/src/graph/traversal/walker.{hpp,cpp} metagraph/tests/graph/traversal/test_walker.cpp metagraph/tests/graph/traversal/test_trie_reference.hpp metagraph/tests/graph/traversal/test_trie_cases.cpp (cases) Co-authored-by: Gunnar Rätsch Co-Authored-By: Claude Opus 5.5 --- metagraph/src/graph/traversal/walker.cpp | 1132 +++++++++++++++++ metagraph/src/graph/traversal/walker.hpp | 246 ++++ .../graph/traversal/test_trie_reference.hpp | 318 +++++ .../tests/graph/traversal/test_walker.cpp | 1049 +++++++++++++++ .../graph/traversal/test_walker_reference.cpp | 339 +++++ .../graph/traversal/walker_test_helpers.hpp | 132 ++ 6 files changed, 3216 insertions(+) create mode 100644 metagraph/src/graph/traversal/walker.cpp create mode 100644 metagraph/src/graph/traversal/walker.hpp create mode 100644 metagraph/tests/graph/traversal/test_trie_reference.hpp create mode 100644 metagraph/tests/graph/traversal/test_walker.cpp create mode 100644 metagraph/tests/graph/traversal/test_walker_reference.cpp create mode 100644 metagraph/tests/graph/traversal/walker_test_helpers.hpp diff --git a/metagraph/src/graph/traversal/walker.cpp b/metagraph/src/graph/traversal/walker.cpp new file mode 100644 index 0000000000..e9afb7eda4 --- /dev/null +++ b/metagraph/src/graph/traversal/walker.cpp @@ -0,0 +1,1132 @@ +#include "walker.hpp" + +#include +#include + +#include +#include + +#include "resolve.hpp" +#include "graph/annotated_dbg.hpp" +#include "graph/representation/canonical_dbg.hpp" +#include "graph/representation/succinct/boss.hpp" +#include "common/seq_tools/reverse_complement.hpp" +#include "common/unix_tools.hpp" +#include "common/vector.hpp" + + +namespace mtg { +namespace graph { +namespace traversal { + +const char* to_string(SeedDerivationError::Cause cause) { + switch (cause) { + case SeedDerivationError::NO_CARRIER: return "no_carrier"; + case SeedDerivationError::TOO_WIDE: return "too_wide"; + case SeedDerivationError::TIME_BUDGET: return "time_budget"; + case SeedDerivationError::AMBIGUOUS_HEADER: return "ambiguous_header"; + } + return "unknown"; +} + +namespace { + +constexpr char kSentinel = boss::BOSS::kSentinel; +// the k-mers whose rows a derivation decodes at once; the narrowest of the first window +// starts the intersection +constexpr size_t kDeriveWindow = 64; + +// One live label on one path head +struct Entry { + LabelId label = 0; + uint32_t branches = 0; + EndReason reason = EndReason::DEAD_END; // set when the label ends +}; +// sorted by label id +using State = std::vector; + +std::vector labels_of(const State &state) { + std::vector out; + out.reserve(state.size()); + for (const Entry &e : state) { + out.push_back(e.label); + } + return out; +} + +// A live path head +struct Item { + size_t segment = 0; + node_index node = npos; + std::string kmer; // spelling of |node| (uppercase, natural orientation) + uint64_t ext_bp = 0; + State state; + uint32_t splits = 0; + size_t path_id = 0; +}; + +struct Succ { + node_index node; + char ch; +}; + +// An edge ((k+1)-mer): 2-bit packed when it fits 64 bits, otherwise a 128-bit FNV-1a pair +// (two distinct edges colliding on it are treated as one: conservative, it blocks) +struct EdgeKey { + uint64_t lo; + uint64_t hi; + bool operator==(const EdgeKey &o) const { return lo == o.lo && hi == o.hi; } +}; +struct EdgeKeyHash { + size_t operator()(const EdgeKey &k) const { + return k.lo ^ (k.hi * 0x9E3779B97F4A7C15ULL) ^ (k.hi >> 29); + } +}; +// one use of an edge: the segment that took it and the orientation bit +struct EdgeUse { + size_t segment; + bool rc; +}; +// the same use keyed the other way round, for probing "did THIS segment take the edge" +struct EdgeSegKey { + EdgeKey edge; + size_t segment; + bool operator==(const EdgeSegKey &o) const { return edge == o.edge && segment == o.segment; } +}; +struct EdgeSegKeyHash { + size_t operator()(const EdgeSegKey &k) const { + return EdgeKeyHash()(k.edge) ^ (k.segment * 0xC2B2AE3D27D4EB4FULL); + } +}; + +// A successor under evaluation +struct Cand { + const Succ *succ = nullptr; + std::string kmer; + State state; // the head's labels that carry the successor + bool skipped = false; // a hairpin, not followed + bool blocked = false; + EndReason block_reason = EndReason::DEAD_END; + EdgeKey edge { 0, 0 }; + bool rc = false; + bool followed = false; + bool split_limited = false; // refused by the path's split allowance + bool admissible() const { return !skipped && !blocked; } +}; + +// a blocked successor's reason, the strongest of several winning +uint8_t block_rank(EndReason reason) { + switch (reason) { + case EndReason::REACHED_SEED: return 3; + case EndReason::EDGE_REUSE_RC: return 2; + case EndReason::EDGE_REUSE: return 1; + default: return 0; + } +} + +struct LeafInfo { + bool is_leaf = false; + std::optional path_reason; + std::vector end_labels; +}; + +struct ArmState { + Arm arm = Arm::RIGHT; + ArmResult result; + std::vector frontier; + std::vector next; + std::vector walk_seq; // per segment, in walking order + std::vector leaves; // per segment + // The per-path edge-reuse record, kept both ways: every use of an edge (scanned when + // the edge has few uses) and, per (edge, segment), the orientation bits 1 << rc of that + // segment's uses (probed per ancestor when it has many: every live path through a + // shared region records the same edge, and scanning those uses once per path would be + // quadratic in the number of paths) + tsl::hopscotch_map, EdgeKeyHash> used_edges; + tsl::hopscotch_map used_by_segment; + // the ancestors of |marked_segment|, as marks and as a list + std::vector visit_mark; + uint32_t visit_epoch = 0; + size_t marked_segment = SIZE_MAX; + std::vector ancestors; + size_t finished_leaves = 0; + size_t next_path_id = 1; + bool stopped = false; + // extension depth of the first head not expanded: its level is partial + uint64_t boundary = std::numeric_limits::max(); +}; + + +class Walker { + public: + Walker(const LabelOracle &oracle, const Seed &seed, const Strategy &strategy) + : oracle_(oracle), graph_(oracle.graph()), seed_(seed), strategy_(strategy), + k_(graph_.get_k()), regime_(oracle.regime()), + canonical_(dynamic_cast(&graph_)) { + size_t n = 0; + code_.fill(255); + for (char c : graph_.alphabet()) { + if (c != kSentinel) + code_[static_cast(c)] = n++; + } + bits_ = 1; + while ((size_t(1) << bits_) < n) { + ++bits_; + } + packable_ = (k_ + 1) * bits_ <= 64; + } + + SeedResult run(); + + private: + const LabelOracle &oracle_; + const DeBruijnGraph &graph_; + const Seed &seed_; + const Strategy &strategy_; + const size_t k_; + const Regime regime_; + const CanonicalDBG *canonical_; + Timer timer_; + + std::array code_; + size_t bits_; + bool packable_; + + SeedResult result_; + std::string seed_upper_; + std::vector nodes_; + tsl::hopscotch_set seed_nodes_; + std::unique_ptr query_; + std::array arms_; + uint64_t steps_total_ = 0; + bool seed_stopped_ = false; + double cap_demand_ = 0; + + double elapsed_ms() const { return timer_.elapsed() * 1000; } + bool time_exceeded() const { + return strategy_.time_budget_ms <= 0 || elapsed_ms() >= strategy_.time_budget_ms; + } + + // seed phase + void validate_seed(); + void derive_seed_labels(const std::vector &keys, std::vector *refs); + void init_arm(ArmState &arm); + + // graph access + std::vector successors(Arm arm, node_index node, const std::string &kmer) const; + void succ_kmer(Arm arm, const std::string &kmer, char c, std::string *out) const; + void step_kmer(Arm arm, const std::string &kmer, char c, std::string *out) const; + bool is_hairpin(const std::string &kmer_u, const std::string &kmer_v, + const std::string &step) const; + std::pair edge_key(const std::string &step) const; + + // bookkeeping + size_t new_segment(ArmState &arm, std::vector parents, uint64_t from_bp, + std::vector labels_start); + void mark_ancestors(ArmState &arm, size_t segment); + bool is_ancestor_or_self(const ArmState &arm, size_t anc, size_t segment) const; + void record_edge(ArmState &arm, size_t segment, const Cand &c); + void end_label(ArmState &arm, Item &item, Entry &e, EndReason reason, const char *text); + void finish_path(ArmState &arm, Item &item, std::optional path_reason); + void censor_item(ArmState &arm, Item &item, EndReason reason); + + // the walk + void run_level(ArmState &arm, uint64_t depth); + void check_structure(ArmState &arm, const Item &item, Cand &c); + std::optional cap_check(const ArmState &arm, size_t nf, size_t remaining); + std::optional process_item(ArmState &arm, Item &item, + const std::vector &succs, + const LabelQuery::NodeHits *hits, + size_t remaining_in_level); + void stop_arm(ArmState &arm, EndReason reason, std::vector *items, size_t from); + void stop_frontier(ArmState &arm, EndReason reason); + void trip(ArmState &arm, EndReason reason, std::vector &items, size_t index); + void finalize(ArmState &arm); +}; + + +/******************************** seed phase ********************************/ + +void Walker::validate_seed() { + const std::string &seq = seed_.sequence; + if (seq.size() < k_) + throw std::invalid_argument("Seed shorter than k"); + + seed_upper_ = seq; +#if ! _DNA_CASE_SENSITIVE_GRAPH + for (char &c : seed_upper_) { + c = std::toupper(static_cast(c)); + } +#endif + const std::string &alphabet = graph_.alphabet(); + for (size_t i = 0; i < seed_upper_.size(); ++i) { + char c = seed_upper_[i]; + if (c == kSentinel || alphabet.find(c) == std::string::npos) { + throw std::invalid_argument("Invalid character '" + std::string(1, seq[i]) + + "' at position " + std::to_string(i) + " of the seed"); + } + } + + nodes_ = map_to_nodes_sequentially(graph_, seed_upper_); + assert(nodes_.size() == seq.size() - k_ + 1); + if (std::find(nodes_.begin(), nodes_.end(), npos) != nodes_.end()) { + const auto runs = runs_of(nodes_.size(), [&](uint64_t i) { return nodes_[i] != npos; }); + uint64_t present = 0; + for (const KmerInterval &run : runs) { + present += run.size(); + } + throw SeedNotInGraph("Seed is not fully present in the graph; graph runs: " + + encode_runs(runs, nodes_.size()), nodes_.size(), present); + } + result_.seed_id = seed_.seed_id; + result_.length_bp = seq.size(); + result_.num_kmers = nodes_.size(); + + // the seed's nodes, in both orientations in the canonical regimes: a step into any of + // them re-enters the seed + for (node_index n : nodes_) { + seed_nodes_.insert(n); + } + if (regime_ != Regime::BASIC) { + std::string rc = seed_upper_; + ::reverse_complement(rc); + for (node_index n : map_to_nodes_sequentially(graph_, rc)) { + if (n != npos) + seed_nodes_.insert(n); + } + } + + const std::vector keys = oracle_.keys_of_sequence(seed_upper_); + std::vector refs; + if (seed_.labels.empty()) { + // every derived label carries every seed k-mer + derive_seed_labels(keys, &refs); + result_.label_dict = refs; + } else { + tsl::hopscotch_set named; + for (const auto &name : seed_.labels) { + if (!named.insert(name).second) + throw std::invalid_argument("Duplicate seed label '" + name + "'"); + refs.push_back(oracle_.resolve_label(name)); + } + // a label that does not carry every seed k-mer is dropped, with the runs it carries + const auto hits = LabelQuery(oracle_, refs).fetch(keys); + std::vector>> runs(refs.size()); + for (uint64_t i = 0; i < hits.size(); ++i) { + for (LabelId l : hits[i]) { + if (runs[l].empty() || runs[l].back().second != i) { + runs[l].emplace_back(i, i + 1); + } else { + runs[l].back().second = i + 1; + } + } + } + for (LabelId l = 0; l < refs.size(); ++l) { + if (runs[l].size() == 1 && runs[l][0] == std::make_pair(uint64_t(0), keys.size())) { + result_.label_dict.push_back(refs[l]); + } else { + result_.dropped_labels.push_back({ refs[l].name, std::move(runs[l]) }); + } + } + if (result_.label_dict.empty()) + throw std::invalid_argument("No seed label supports every k-mer of the seed"); + result_.labels_supporting_total = result_.label_dict.size(); + } + + std::vector names; + for (const LabelRef &ref : result_.label_dict) { + names.push_back(ref.name); + } + result_.validated_seed_id = make_seed_id(seed_upper_, regime_ != Regime::BASIC, names); + result_.seed_id_mismatch = !seed_.seed_id.empty() + && seed_.seed_id != result_.validated_seed_id; +} + +void Walker::derive_seed_labels(const std::vector &keys, + std::vector *refs) { + using Key = std::pair; + + const LabelKind kind = strategy_.seed_label_kind.value_or( + oracle_.coord_to_header() ? LabelKind::HEADER : LabelKind::COLUMN); + if (kind == LabelKind::HEADER && !oracle_.coord_to_header()) { + throw std::invalid_argument("Deriving sequence header labels from the seed requires a " + "CoordToHeader index; set seed_label_kind to \"column\""); + } + // a seed of exactly k bases has no intersection to narrow the candidates: bound what + // the first row may materialise + const size_t max_candidates = std::max( + size_t(1) << 16, 64 * std::min(strategy_.max_seed_labels, size_t(1) << 40)); + + // the labels carrying every k-mer read so far, ascending + std::vector live, next, cur; + size_t read = 0; + for (size_t begin = 0; begin < keys.size(); begin += kDeriveWindow) { + const size_t end = std::min(keys.size(), begin + kDeriveWindow); + std::vector distinct; + for (size_t i = begin; i < end; ++i) { + distinct.push_back(AnnotatedDBG::graph_to_anno_index(keys[i])); + } + std::sort(distinct.begin(), distinct.end()); + distinct.erase(std::unique(distinct.begin(), distinct.end()), distinct.end()); + std::vector plain; + std::vector tuples; + if (kind == LabelKind::HEADER) { + tuples = oracle_.get_row_tuples(distinct); + } else { + plain = oracle_.get_rows(distinct); + } + auto row_of = [&](size_t i) { + return static_cast( + std::lower_bound(distinct.begin(), distinct.end(), + AnnotatedDBG::graph_to_anno_index(keys[i])) + - distinct.begin()); + }; + // the size of a row: its columns, and for a tuple row its coordinates + auto cost_of = [&](size_t i) { + const size_t r = row_of(i); + if (kind == LabelKind::COLUMN) + return plain[r].size(); + size_t n = 0; + for (const auto &entry : tuples[r]) { + n += 1 + entry.second.size(); + } + return n; + }; + + std::vector order; + for (size_t i = begin; i < end; ++i) { + order.push_back(i); + } + if (!read) { + // the narrowest k-mer of the first window starts the intersection + std::vector cost; + for (size_t i : order) { + cost.push_back(cost_of(i)); + } + const size_t best = std::min_element(cost.begin(), cost.end()) - cost.begin(); + std::iter_swap(order.begin(), order.begin() + best); + if (cost[best] > max_candidates) { + throw SeedDerivationError(SeedDerivationError::TOO_WIDE, + "The permitted set cannot be derived from this seed: the narrowest of " + "its first " + std::to_string(end - begin) + " k-mers alone has " + + std::to_string(cost[best]) + " annotation entries, over the limit of " + + std::to_string(max_candidates) + ". Start the seed in a less " + "repetitive k-mer, or name the labels explicitly.", + static_cast(max_candidates), static_cast(cost[best])); + } + } + + for (size_t i : order) { + const size_t r = row_of(i); + cur.clear(); + if (kind == LabelKind::COLUMN) { + for (Column c : plain[r]) { + cur.emplace_back(c, 0); + } + } else { + for (const auto &[c, coords] : tuples[r]) { + // only the columns still live can contribute + if (read && !std::binary_search(live.begin(), live.end(), Key(c, 0), + [](const Key &a, const Key &b) { + return a.first < b.first; + })) { + continue; + } + for (uint64_t coord : coords) { + cur.emplace_back(c, oracle_.map_coord(c, coord).first); + } + } + } + std::sort(cur.begin(), cur.end()); + cur.erase(std::unique(cur.begin(), cur.end()), cur.end()); + if (!read) { + live.swap(cur); + } else { + next.clear(); + std::set_intersection(live.begin(), live.end(), cur.begin(), cur.end(), + std::back_inserter(next)); + live.swap(next); + } + ++read; + if (live.empty()) { + throw SeedDerivationError(SeedDerivationError::NO_CARRIER, + "No label supports every k-mer of the seed (the " + "permitted set derived from the seed is empty at " + "k-mer " + std::to_string(i) + " of " + + std::to_string(keys.size()) + ")", + static_cast(keys.size()), + static_cast(read)); + } + if (read < keys.size() && strategy_.time_budget_ms > 0 + && elapsed_ms() >= strategy_.time_budget_ms) { + throw SeedDerivationError(SeedDerivationError::TIME_BUDGET, + "The time budget (" + std::to_string(strategy_.time_budget_ms) + + " ms) ran out deriving the permitted set from the seed, after " + + std::to_string(read) + " of its " + std::to_string(keys.size()) + + " k-mers; raise bounds.time_budget_ms or name the labels explicitly", + strategy_.time_budget_ms, static_cast(read)); + } + } + } + + result_.labels_from_seed = true; + result_.labels_supporting_total = live.size(); + auto name_of = [&](const Key &key) -> const std::string& { + return kind == LabelKind::COLUMN ? oracle_.column_name(key.first) + : oracle_.header_name(key.first, key.second); + }; + if (kind == LabelKind::HEADER) { + // A derived list must be resubmittable verbatim as an explicit one, which resolves a + // name to a column first, then to the first column holding the header + for (const Key &key : live) { + const std::string &name = name_of(key); + const LabelRef back = oracle_.resolve_label(name); // the header exists: no throw + if (back.kind != LabelKind::HEADER || back.column != key.first + || back.seq_id != key.second) { + throw SeedDerivationError(SeedDerivationError::AMBIGUOUS_HEADER, + "The labels derived from the seed are not resubmittable: the sequence " + "header '" + name + "' " + + (back.kind == LabelKind::COLUMN + ? "is also the name of an annotation column" + : "also occurs in another annotation column") + + ". Set seed_label_kind to \"column\" or name the labels explicitly.", + 0, 0, name); + } + } + } + if (live.size() > strategy_.max_seed_labels) { + uint64_t digest = 0xcbf29ce484222325ULL; + for (size_t i = strategy_.max_seed_labels; i < live.size(); ++i) { + digest = fnv1a64(name_of(live[i]) + "\n", digest); + } + result_.labels_dropped = live.size() - strategy_.max_seed_labels; + result_.labels_dropped_digest = hex64(digest); + live.resize(strategy_.max_seed_labels); + } + for (const Key &key : live) { + LabelRef ref; + ref.kind = kind; + ref.column = key.first; + ref.seq_id = kind == LabelKind::HEADER ? key.second : 0; + ref.name = name_of(key); + refs->push_back(std::move(ref)); + } +} + +void Walker::init_arm(ArmState &arm) { + arm.result.arm = arm.arm; + arm.result.requested = strategy_.direction == Strategy::BOTH + || (strategy_.direction == Strategy::LEFT) == (arm.arm == Arm::LEFT); + if (!arm.result.requested) + return; + + Item root; + root.node = arm.arm == Arm::RIGHT ? nodes_.back() : nodes_.front(); + root.kmer = arm.arm == Arm::RIGHT ? seed_upper_.substr(seed_upper_.size() - k_) + : seed_upper_.substr(0, k_); + for (LabelId l = 0; l < result_.label_dict.size(); ++l) { + root.state.push_back(Entry{ l, 0 }); + } + root.segment = new_segment(arm, {}, 0, labels_of(root.state)); + arm.frontier.push_back(std::move(root)); +} + + +/******************************* graph access *******************************/ + +std::vector Walker::successors(Arm arm, node_index node, const std::string &kmer) const { + std::vector out; + auto cb = [&](node_index n, char c) { + if (c != kSentinel && n != npos) + out.push_back({ n, c }); + }; + try { + if (arm == Arm::RIGHT) { + if (canonical_) { + canonical_->call_outgoing_kmers(node, kmer, cb); + } else { + graph_.call_outgoing_kmers(node, cb); + } + } else { + if (canonical_) { + canonical_->call_incoming_kmers(node, kmer, cb); + } else { + graph_.call_incoming_kmers(node, cb); + } + } + } catch (const std::runtime_error &e) { + // CanonicalDBG throws an empty runtime_error on an inconsistent primary graph + throw std::runtime_error("Inconsistent primary graph at node " + std::to_string(node) + + " (" + kmer + ") on arm " + to_string(arm) + + (e.what()[0] ? std::string(": ") + e.what() : std::string())); + } + std::sort(out.begin(), out.end(), [](const Succ &a, const Succ &b) { return a.ch < b.ch; }); + return out; +} + +// the k-mer of the successor reached by |c| (|out| must not alias |kmer|) +void Walker::succ_kmer(Arm arm, const std::string &kmer, char c, std::string *out) const { + out->clear(); + if (arm == Arm::RIGHT) { + out->append(kmer, 1, k_ - 1); + out->push_back(c); + } else { + out->push_back(c); + out->append(kmer, 0, k_ - 1); + } +} + +// the (k+1)-mer of a step in natural orientation +void Walker::step_kmer(Arm arm, const std::string &kmer, char c, std::string *out) const { + out->clear(); + if (arm == Arm::RIGHT) { + out->append(kmer); + out->push_back(c); + } else { + out->push_back(c); + out->append(kmer); + } +} + +// A self-reverse-complementary step: the (k+1)-mer is its own reverse complement, or for +// even k either k-mer is +bool Walker::is_hairpin(const std::string &kmer_u, const std::string &kmer_v, + const std::string &step) const { + auto palindromic = [](const std::string &s) { + std::string rc = s; + ::reverse_complement(rc); + return rc == s; + }; + return palindromic(step) || (k_ % 2 == 0 && (palindromic(kmer_v) || palindromic(kmer_u))); +} + +// the edge of a step, canonical in the canonical regimes, and whether it was flipped +std::pair Walker::edge_key(const std::string &step) const { + std::string rc; + std::string_view s = step; + bool flip = false; + if (regime_ != Regime::BASIC) { + rc = step; + ::reverse_complement(rc); + if (rc < step) { + s = rc; + flip = true; + } + } + EdgeKey key { 0, 0 }; + if (packable_) { + for (char c : s) { + key.lo = (key.lo << bits_) | code_[static_cast(c)]; + } + } else { + key.lo = fnv1a64(s); + key.hi = fnv1a64(s, 0x84222325cbf29ce4ULL); + } + return { key, flip }; +} + + +/******************************** bookkeeping *******************************/ + +size_t Walker::new_segment(ArmState &arm, std::vector parents, uint64_t from_bp, + std::vector labels_start) { + Segment seg; + seg.id = arm.result.segments.size(); + seg.parents = std::move(parents); + seg.from_bp = from_bp; + seg.labels_start = std::move(labels_start); + arm.result.segments.push_back(std::move(seg)); + arm.walk_seq.emplace_back(); + arm.leaves.emplace_back(); + return arm.result.segments.size() - 1; +} + +// Mark the ancestors of |segment| (its chain of parents). The parents of a segment never +// change, so the marks stay valid until another segment is marked. +void Walker::mark_ancestors(ArmState &arm, size_t segment) { + if (arm.marked_segment == segment) + return; + const auto &segments = arm.result.segments; + if (arm.visit_mark.size() < segments.size()) + arm.visit_mark.resize(segments.size(), 0); + ++arm.visit_epoch; + arm.ancestors.clear(); + for (size_t s = segment; !segments[s].parents.empty(); ) { + s = segments[s].parents[0]; + arm.visit_mark[s] = arm.visit_epoch; + arm.ancestors.push_back(s); + } + arm.marked_segment = segment; +} + +// requires mark_ancestors(arm, segment) +bool Walker::is_ancestor_or_self(const ArmState &arm, size_t anc, size_t segment) const { + assert(arm.marked_segment == segment); + // segments are numbered in creation order + return anc == segment || (anc < segment && arm.visit_mark[anc] == arm.visit_epoch); +} + +void Walker::record_edge(ArmState &arm, size_t segment, const Cand &c) { + arm.used_edges[c.edge].push_back(EdgeUse{ segment, c.rc }); + arm.used_by_segment[EdgeSegKey{ c.edge, segment }] |= static_cast(1u << c.rc); +} + +void Walker::end_label(ArmState &arm, Item &item, Entry &e, EndReason reason, + const char *text) { + e.reason = reason; + arm.result.segments[item.segment].label_ends.push_back( + LabelEnd{ item.ext_bp, e.label, reason, text }); +} + +// |item| is a leaf: every label on it has ended +void Walker::finish_path(ArmState &arm, Item &item, std::optional path_reason) { + arm.result.segments[item.segment].labels_end = labels_of(item.state); + LeafInfo &leaf = arm.leaves[item.segment]; + leaf.is_leaf = true; + leaf.path_reason = path_reason; + for (const Entry &e : item.state) { + leaf.end_labels.push_back(PathEnd{ e.label, e.reason, e.branches }); + } + arm.finished_leaves++; +} + +void Walker::censor_item(ArmState &arm, Item &item, EndReason reason) { + for (Entry &e : item.state) { + end_label(arm, item, e, reason, ""); + } + finish_path(arm, item, reason); +} + + +/********************************** the walk ********************************/ + +void Walker::check_structure(ArmState &arm, const Item &item, Cand &c) { + std::string step; + step_kmer(arm.arm, item.kmer, c.succ->ch, &step); + if (regime_ != Regime::BASIC && is_hairpin(item.kmer, c.kmer, step)) { + c.skipped = true; + return; + } + if (seed_nodes_.count(c.succ->node)) { + c.blocked = true; + c.block_reason = EndReason::REACHED_SEED; + return; + } + std::tie(c.edge, c.rc) = edge_key(step); + auto it = arm.used_edges.find(c.edge); + if (it == arm.used_edges.end()) + return; + // Was the edge taken by this path (the head's segment or an ancestor), and in which + // orientation? Either scan the edge's uses or probe (edge, segment) for the path's own + // segments, whichever are fewer: the same uses, so the same answer. + mark_ancestors(arm, item.segment); + const std::vector &uses = it->second; + bool same = false, opposite = false; + if (uses.size() <= arm.ancestors.size() + 1) { + for (const EdgeUse &u : uses) { + if (is_ancestor_or_self(arm, u.segment, item.segment)) + (u.rc == c.rc ? same : opposite) = true; + } + } else { + auto probe = [&](size_t segment) { + auto jt = arm.used_by_segment.find(EdgeSegKey{ c.edge, segment }); + if (jt == arm.used_by_segment.end()) + return; + same |= (jt->second >> c.rc) & 1; + opposite |= (jt->second >> !c.rc) & 1; + }; + probe(item.segment); + for (size_t a : arm.ancestors) { + probe(a); + } + } + if (opposite) { + c.blocked = true; + c.block_reason = EndReason::EDGE_REUSE_RC; + } else if (same) { + c.blocked = true; + c.block_reason = EndReason::EDGE_REUSE; + } +} + +// The caps, decided before anything of a head is committed: |nf| successors would be +// followed, |remaining| heads of the level are still to come after this one +std::optional Walker::cap_check(const ArmState &arm, size_t nf, size_t remaining) { + if (!nf) + return std::nullopt; + auto over = [&](uint64_t demand, EndReason reason) { + cap_demand_ = static_cast(demand); + return reason; + }; + if (steps_total_ + nf > strategy_.max_steps) + return over(steps_total_ + nf, EndReason::MAX_STEPS); + if (arm.result.output_bp + nf > strategy_.max_output_bp) + return over(arm.result.output_bp + nf, EndReason::MAX_OUTPUT); + const size_t live_after = remaining + arm.next.size() + nf; + if (nf >= 2 && arm.finished_leaves + live_after > strategy_.max_paths) + return over(arm.finished_leaves + live_after, EndReason::MAX_PATHS); + if (live_after > strategy_.max_live_paths) + return over(live_after, EndReason::MAX_LIVE_PATHS); + return std::nullopt; +} + +std::optional Walker::process_item(ArmState &arm, Item &item, + const std::vector &succs, + const LabelQuery::NodeHits *hits, + size_t remaining_in_level) { + const uint64_t at = item.ext_bp; + if (succs.empty()) { + for (Entry &e : item.state) { + end_label(arm, item, e, EndReason::DEAD_END, ""); + } + finish_path(arm, item, std::nullopt); + return std::nullopt; + } + + // ---- the head's labels on each successor, and the walk rule + std::vector cands(succs.size()); + for (size_t i = 0; i < succs.size(); ++i) { + Cand &c = cands[i]; + c.succ = &succs[i]; + succ_kmer(arm.arm, item.kmer, succs[i].ch, &c.kmer); + // both sorted by label + auto h = hits[i].begin(); + for (const Entry &e : item.state) { + while (h != hits[i].end() && *h < e.label) { + ++h; + } + if (h != hits[i].end() && *h == e.label) + c.state.push_back(e); + } + check_structure(arm, item, c); + } + + // ---- per-lineage branching: a label on several admissible successors continues on + // all of them within its branch allowance, else it ends on all of them (BRANCH) + tsl::hopscotch_map continuations; + for (const Cand &c : cands) { + if (!c.admissible()) + continue; + for (const Entry &e : c.state) { + continuations[e.label]++; + } + } + tsl::hopscotch_set excluded; + for (const Entry &src : item.state) { + auto it = continuations.find(src.label); + if (it != continuations.end() && it->second >= 2 + && (strategy_.max_label_branches == Strategy::kUnlimited + || src.branches + 1 > strategy_.max_label_branches)) { + excluded.insert(src.label); + } + } + bool ambiguous = false; + std::vector followed; + for (Cand &c : cands) { + if (!c.admissible()) + continue; + if (!excluded.empty()) { + c.state.erase(std::remove_if(c.state.begin(), c.state.end(), [&](const Entry &e) { + return excluded.count(e.label); + }), c.state.end()); + } + for (Entry &e : c.state) { + if (continuations[e.label] >= 2) { + e.branches++; + ambiguous = true; + } + } + if (!c.state.empty()) { + c.followed = true; + followed.push_back(&c); + } + } + if (followed.size() >= 2 && item.splits + 1 > strategy_.max_splits_per_path) { + for (Cand *c : followed) { + c->followed = false; + c->split_limited = true; + } + followed.clear(); + } + const size_t nf = followed.size(); + + // ---- the caps, before anything is committed + if (auto cap = cap_check(arm, nf, remaining_in_level)) + return cap; + + // ---- the labels that end here, in state order + for (Entry &src : item.state) { + const LabelId l = src.label; + bool continues = false, split_limited = false, hairpin = false; + uint8_t blocked = 0; + for (const Cand &c : cands) { + const bool has = std::binary_search(c.state.begin(), c.state.end(), src, + [](const Entry &a, const Entry &b) { return a.label < b.label; }); + if (!has) + continue; + continues |= c.followed; + split_limited |= c.split_limited; + hairpin |= c.skipped; + if (c.blocked) + blocked = std::max(blocked, block_rank(c.block_reason)); + } + if (continues) + continue; + if (excluded.count(l)) { + end_label(arm, item, src, EndReason::BRANCH, ""); + } else if (split_limited) { + end_label(arm, item, src, EndReason::BRANCH, "split_limit"); + } else if (blocked) { + end_label(arm, item, src, blocked == 3 ? EndReason::REACHED_SEED + : blocked == 2 ? EndReason::EDGE_REUSE_RC + : EndReason::EDGE_REUSE, ""); + } else if (hairpin) { + end_label(arm, item, src, EndReason::DEAD_END, "hairpin"); + } else { + end_label(arm, item, src, EndReason::LABEL_LOST, ""); + } + } + + if (!nf) { + finish_path(arm, item, std::nullopt); + return std::nullopt; + } + + // ---- steps + arm.result.steps += nf; + arm.result.output_bp += nf; + steps_total_ += nf; + + if (nf == 1) { + Cand &c = *followed[0]; + arm.walk_seq[item.segment].push_back(c.succ->ch); + arm.result.segments[item.segment].length_bp++; + record_edge(arm, item.segment, c); + item.node = c.succ->node; + item.kmer = std::move(c.kmer); + item.ext_bp = at + 1; + item.state = std::move(c.state); + arm.next.push_back(std::move(item)); + return std::nullopt; + } + + // ---- split + Split split { at, item.segment, {}, ambiguous }; + arm.result.segments[item.segment].labels_end = labels_of(item.state); + for (Cand *c : followed) { + size_t child = new_segment(arm, { item.segment }, at, labels_of(c->state)); + arm.result.segments[item.segment].children.push_back(child); + split.children.push_back(child); + arm.walk_seq[child].push_back(c->succ->ch); + arm.result.segments[child].length_bp = 1; + record_edge(arm, child, *c); + Item ni; + ni.segment = child; + ni.node = c->succ->node; + ni.kmer = std::move(c->kmer); + ni.ext_bp = at + 1; + ni.state = std::move(c->state); + ni.splits = item.splits + 1; + ni.path_id = arm.next_path_id++; + arm.next.push_back(std::move(ni)); + } + arm.result.splits.push_back(std::move(split)); + return std::nullopt; +} + +// Stop the arm at a cap: the heads from items[from] on and every head of the next level are +// censored; a head that already reached the extension length ends with MAX_EXTENSION +void Walker::stop_arm(ArmState &arm, EndReason reason, std::vector *items, size_t from) { + const size_t heads = arm.next.size() + (items ? items->size() - from : 0); + if (!heads) + return; + const uint64_t radius = strategy_.max_extension_bp; + size_t live = 0; + for (const Item &item : arm.next) { + live += item.ext_bp < radius; + } + if (items) { + for (size_t j = from; j < items->size(); ++j) { + live += (*items)[j].ext_bp < radius; + } + } + const Item &first = items && from < items->size() ? (*items)[from] : arm.next.front(); + if (!arm.result.cap_trigger) + arm.result.cap_trigger = CapTrigger{ reason, first.ext_bp, cap_demand_ }; + arm.boundary = std::min(arm.boundary, first.ext_bp); + arm.result.frontier_live_paths = live; + arm.result.status = ArmResult::TRUNCATED; + auto censor = [&](Item &item) { + censor_item(arm, item, item.ext_bp >= radius ? EndReason::MAX_EXTENSION : reason); + }; + if (items) { + for (size_t j = from; j < items->size(); ++j) { + censor((*items)[j]); + } + } + for (Item &item : arm.next) { + censor(item); + } + arm.next.clear(); + arm.frontier.clear(); + arm.stopped = true; +} + +void Walker::stop_frontier(ArmState &arm, EndReason reason) { + if (arm.frontier.empty()) + return; + std::vector rest; + for (Item &item : arm.frontier) { + if (item.ext_bp >= strategy_.max_extension_bp) { + censor_item(arm, item, EndReason::MAX_EXTENSION); + } else { + rest.push_back(std::move(item)); + } + } + arm.frontier.swap(rest); + if (!arm.frontier.empty()) + stop_arm(arm, reason, &arm.frontier, 0); +} + +// A cap at items[index]: the steps and the time stop the whole seed, the others the arm +void Walker::trip(ArmState &arm, EndReason reason, std::vector &items, size_t index) { + stop_arm(arm, reason, &items, index); + if (reason == EndReason::MAX_STEPS || reason == EndReason::TIME_BUDGET) { + seed_stopped_ = true; + for (ArmState &other : arms_) { + if (&other != &arm) + stop_frontier(other, reason); + } + } +} + +void Walker::run_level(ArmState &arm, uint64_t depth) { + std::vector items; + items.swap(arm.frontier); + std::sort(items.begin(), items.end(), + [](const Item &a, const Item &b) { return a.path_id < b.path_id; }); + + // the successors of every head below the extension length, and their labels in one fetch + std::vector> succs(items.size()); + std::vector keys; + std::string kmer; + for (size_t i = 0; i < items.size(); ++i) { + if (items[i].ext_bp >= strategy_.max_extension_bp) + continue; + succs[i] = successors(arm.arm, items[i].node, items[i].kmer); + for (const Succ &s : succs[i]) { + // the spelling is only needed by native canonical graphs + if (regime_ == Regime::CANONICAL) + succ_kmer(arm.arm, items[i].kmer, s.ch, &kmer); + keys.push_back(oracle_.key_of(s.node, kmer)); + } + } + const std::vector hits = query_->fetch(keys); + + size_t offset = 0; + for (size_t i = 0; i < items.size(); ++i) { + Item &item = items[i]; + if (item.ext_bp >= strategy_.max_extension_bp) { + censor_item(arm, item, EndReason::MAX_EXTENSION); + continue; + } + // the deadline is read before every head; never at depth 0, so that a time budget + // of 0 still delivers the result complete to the first level + if (depth > 0 && time_exceeded()) { + cap_demand_ = elapsed_ms(); + trip(arm, EndReason::TIME_BUDGET, items, i); + return; + } + auto cap = process_item(arm, item, succs[i], hits.data() + offset, items.size() - i - 1); + offset += succs[i].size(); + if (cap) { + trip(arm, *cap, items, i); + return; + } + } + arm.frontier.swap(arm.next); +} + +void Walker::finalize(ArmState &arm) { + ArmResult &res = arm.result; + // every walk up to the first unexpanded head's depth is present + res.complete_to_bp = std::min(arm.boundary, strategy_.max_extension_bp); + for (size_t s = 0; s < res.segments.size(); ++s) { + Segment &seg = res.segments[s]; + if (strategy_.sequences) { + seg.sequence = std::move(arm.walk_seq[s]); + if (arm.arm == Arm::LEFT) + std::reverse(seg.sequence.begin(), seg.sequence.end()); + } + std::stable_sort(seg.label_ends.begin(), seg.label_ends.end(), + [](const LabelEnd &a, const LabelEnd &b) { return a.at_bp < b.at_bp; }); + } + // one path per leaf, in segment order + for (size_t s = 0; s < res.segments.size(); ++s) { + LeafInfo &leaf = arm.leaves[s]; + if (!leaf.is_leaf) + continue; + PathResult path; + path.id = res.paths.size(); + path.leaf = s; + path.length_bp = res.segments[s].from_bp + res.segments[s].length_bp; + path.path_reason = leaf.path_reason; + path.end_labels = std::move(leaf.end_labels); + res.paths.push_back(std::move(path)); + } +} + +SeedResult Walker::run() { + if (!strategy_.max_seed_labels) + throw std::invalid_argument("max_seed_labels must be positive"); + + validate_seed(); + query_ = std::make_unique(oracle_, result_.label_dict); + + arms_[0].arm = Arm::LEFT; + arms_[1].arm = Arm::RIGHT; + for (ArmState &arm : arms_) { + init_arm(arm); + } + + for (uint64_t depth = 0; !seed_stopped_; ++depth) { + if (std::all_of(arms_.begin(), arms_.end(), + [](const ArmState &arm) { return arm.frontier.empty(); })) + break; + if (depth > 0 && time_exceeded()) { + // only heads below the extension length are censored by the deadline: when every + // remaining head reached it, the walk is complete and nothing was stopped + cap_demand_ = elapsed_ms(); + for (ArmState &arm : arms_) { + stop_frontier(arm, EndReason::TIME_BUDGET); + } + break; + } + for (ArmState &arm : arms_) { + if (!arm.frontier.empty() && !arm.stopped) + run_level(arm, depth); + if (seed_stopped_) + break; + } + } + + for (ArmState &arm : arms_) { + finalize(arm); + } + result_.arms[static_cast(Arm::LEFT)] = std::move(arms_[0].result); + result_.arms[static_cast(Arm::RIGHT)] = std::move(arms_[1].result); + result_.annotation_counters = oracle_.counters(); + result_.elapsed_seconds = elapsed_ms() / 1000; + return std::move(result_); +} + +} // namespace + + +SeedResult traverse_seed(const LabelOracle &oracle, const Seed &seed, const Strategy &strategy) { + return Walker(oracle, seed, strategy).run(); +} + +} // namespace traversal +} // namespace graph +} // namespace mtg diff --git a/metagraph/src/graph/traversal/walker.hpp b/metagraph/src/graph/traversal/walker.hpp new file mode 100644 index 0000000000..d248b9a551 --- /dev/null +++ b/metagraph/src/graph/traversal/walker.hpp @@ -0,0 +1,246 @@ +#ifndef __TRAVERSAL_WALKER_HPP__ +#define __TRAVERSAL_WALKER_HPP__ + +#include +#include +#include +#include +#include +#include + +#include "traversal_types.hpp" +#include "label_oracle.hpp" + + +namespace mtg { +namespace graph { +namespace traversal { + +/** + * Label-constrained extension of a seed. + * + * A seed is a sequence whose k-mers are all in the graph, and a set of permitted labels that + * carry every one of them: named explicitly, or derived as the labels on every seed k-mer. + * From the seed's last k-mer (right arm) and its first (left arm) the walker extends level + * by level, one base per level, keeping per path the labels that carry every k-mer walked + * so far. A label leaves a path when no admissible successor carries it (LABEL_LOST, or + * DEAD_END without successors), when its continuation is blocked by the walk rule + * (EDGE_REUSE, EDGE_REUSE_RC, REACHED_SEED), or when it is on several successors beyond + * its branch allowance (BRANCH). A path splits where its labels continue on several + * successors. The result is a tree of segments per arm. + * + * The walk rule, which makes the set of walks finite on a graph with cycles: a path never + * takes an edge ((k+1)-mer) twice, in either orientation in canonical graphs, never steps + * back into the seed, and never takes a self-reverse-complementary step (a hairpin). + * + * Completeness: levels are synchronous, so a cap trips between two heads of one level. + * ArmResult::complete_to_bp is the extension up to which every admissible walk is present. + */ +struct Strategy { + enum Direction { BOTH, LEFT, RIGHT }; + + // "no limit" for max_label_branches + static constexpr size_t kUnlimited = std::numeric_limits::max(); + + Direction direction = BOTH; + + // Cap on a permitted set DERIVED from the seed: the labels above it (in ascending + // (column, seq_id) order) are reported as SeedResult::labels_dropped and its digest. + size_t max_seed_labels = 1000; + // What a derived label denotes. Unset: HEADER when the index has a CoordToHeader (the + // labels are the indexed sequences), COLUMN otherwise. + std::optional seed_label_kind; + + // How many times a label may continue on more than one successor (0: it ends with + // BRANCH at its first fork; kUnlimited: every label-consistent walk is followed) + size_t max_label_branches = 0; + + uint64_t max_extension_bp = 5000; + double time_budget_ms = 30'000; + + // output + bool sequences = true; + + // Safety caps. Not request fields: they keep a request within bounded memory. + size_t max_splits_per_path = 64; + uint64_t max_steps = 1'000'000; + size_t max_live_paths = 1000; + size_t max_paths = 10'000; + uint64_t max_output_bp = 2'000'000; +}; + +struct Seed { + std::string seed_id; + std::string sequence; + // The permitted labels. Empty: derive them from the seed (the labels on every k-mer). + std::vector labels; +}; + +// A seed with a k-mer the graph does not have +class SeedNotInGraph : public std::invalid_argument { + public: + SeedNotInGraph(const std::string &what, uint64_t kmers, uint64_t present) + : std::invalid_argument(what), kmers_(kmers), present_(present) {} + uint64_t kmers() const { return kmers_; } + uint64_t present() const { return present_; } + + private: + uint64_t kmers_; + uint64_t present_; +}; + +/** + * The permitted set could not be derived from a seed (Seed::labels empty). A property of + * the seed and the data, not of the request: a multi-seed request reports it per seed. + */ +class SeedDerivationError : public std::invalid_argument { + public: + enum Cause { + // no label carries every k-mer of the seed: |limit| the seed's k-mers, |observed| + // the k-mers read when the intersection became empty + NO_CARRIER, + // the narrowest of the first 64 k-mers has more annotation entries than the + // derivation materialises: |limit| that bound, |observed| its entries + TOO_WIDE, + // the time budget ran out before every k-mer was read: |limit| the budget, + // |observed| the k-mers read + TIME_BUDGET, + // a derived header name is not resubmittable as an explicit label: |subject| + AMBIGUOUS_HEADER, + }; + + SeedDerivationError(Cause cause, const std::string &what, double limit = 0, + double observed = 0, std::string subject = "") + : std::invalid_argument(what), cause_(cause), limit_(limit), observed_(observed), + subject_(std::move(subject)) {} + + Cause cause() const { return cause_; } + double limit() const { return limit_; } + double observed() const { return observed_; } + const std::string& subject() const { return subject_; } + + private: + Cause cause_; + double limit_; + double observed_; + std::string subject_; +}; +const char* to_string(SeedDerivationError::Cause cause); + +// A label leaving a path at |at_bp| (bases from the seed boundary) +struct LabelEnd { + uint64_t at_bp = 0; + LabelId label = 0; + EndReason reason = EndReason::DEAD_END; + // why a BRANCH or DEAD_END: "split_limit" (the path's split allowance ran out) or + // "hairpin" (the label's only continuation was a skipped hairpin); otherwise empty + const char *text = ""; +}; + +struct Segment { + size_t id = 0; + std::vector parents; // none for the root, one otherwise + std::vector children; + uint64_t from_bp = 0; + uint64_t length_bp = 0; + std::string sequence; // bases added, natural orientation + // the labels live at the segment's entry node and at its last node + std::vector labels_start; + std::vector labels_end; + std::vector label_ends; // sorted by at_bp +}; + +struct Split { + uint64_t at_bp = 0; // = the shared prefix length from the seed boundary + size_t segment = 0; + std::vector children; + bool ambiguous = false; // some label followed several children +}; + +struct PathEnd { + LabelId label = 0; + EndReason reason = EndReason::DEAD_END; + uint32_t branches = 0; // the label's forks with several continuations taken +}; + +// One leaf of the segment tree; the path is its leaf's chain of parents +struct PathResult { + size_t id = 0; + size_t leaf = 0; // the leaf segment + uint64_t length_bp = 0; + // why the path itself ended when no label decided it: the extension length or a cap + std::optional path_reason; + std::vector end_labels; +}; + +struct CapTrigger { + EndReason reason; + uint64_t at_bp; // depth of the first head not expanded + // what the cap compared at the trip, which exceeded its limit (steps, bases, paths, + // live paths, or elapsed milliseconds) + double demand = 0; +}; + +struct ArmResult { + enum Status { COMPLETE, TRUNCATED }; + + Arm arm = Arm::RIGHT; + bool requested = true; + Status status = COMPLETE; + size_t frontier_live_paths = 0; // heads below the extension length when stopped + std::optional cap_trigger; + // Every admissible walk of at most complete_to_bp bases from the seed boundary is in + // |segments|, and every path that ended before it ended for the reported reason. + // Equals Strategy::max_extension_bp exactly when status == COMPLETE. + uint64_t complete_to_bp = 0; + + std::vector segments; + std::vector splits; + std::vector paths; + + uint64_t steps = 0; + uint64_t output_bp = 0; +}; + +// A seed label that does not carry every seed k-mer +struct DroppedLabel { + std::string name; + std::vector> runs; // the k-mer runs it carries +}; + +struct SeedResult { + std::string seed_id; + std::string validated_seed_id; + bool seed_id_mismatch = false; + uint64_t length_bp = 0; + uint64_t num_kmers = 0; + // the permitted labels; a label id is an index here + std::vector label_dict; + // The labels were derived from the seed. Then |labels_dropped| is how many of them + // max_seed_labels cut, with an FNV-1a-64 hex digest of their names. + // |labels_supporting_total| is how many labels carry every seed k-mer. + bool labels_from_seed = false; + size_t labels_supporting_total = 0; + size_t labels_dropped = 0; + std::string labels_dropped_digest; + std::vector dropped_labels; + std::array arms; // indexed by Arm + LabelOracle::Counters annotation_counters; + double elapsed_seconds = 0; +}; + +/** + * Validate the seed and extend it under |strategy|. + * Throws SeedNotInGraph for a seed with a k-mer the graph does not have, + * SeedDerivationError when an empty |seed.labels| yields no permitted set, and + * std::invalid_argument for invalid seeds, labels or strategies. + */ +SeedResult traverse_seed(const LabelOracle &oracle, + const Seed &seed, + const Strategy &strategy); + +} // namespace traversal +} // namespace graph +} // namespace mtg + +#endif // __TRAVERSAL_WALKER_HPP__ diff --git a/metagraph/tests/graph/traversal/test_trie_reference.hpp b/metagraph/tests/graph/traversal/test_trie_reference.hpp new file mode 100644 index 0000000000..0ba5fd6d94 --- /dev/null +++ b/metagraph/tests/graph/traversal/test_trie_reference.hpp @@ -0,0 +1,318 @@ +#ifndef __TEST_TRIE_REFERENCE_HPP__ +#define __TEST_TRIE_REFERENCE_HPP__ + +#include +#include +#include +#include +#include +#include + +#include + +#include "graph/representation/base/sequence_graph.hpp" +#include "graph/traversal/traversal_types.hpp" +#include "common/seq_tools/reverse_complement.hpp" + +/** + * A string-level model of an annotated index and of the walk rule, built from the records + * a fixture is made of and nothing else: no graph object, no annotation matrix, no walker + * code. label_trie() gives the walks the rule admits over ONE label's k-mers — the label's + * maximal label-consistent walks, which is what the walker claims for that label with an + * unlimited branch allowance. + * + * Orientation: a k-mer and its reverse complement are one node in canonical and primary + * graphs, so there both orientations of every record k-mer are in the model; in a basic + * graph only the record's orientation is. Walks are spelled OUTWARD from the seed boundary + * (the left arm's bases in order of distance from the seed). + */ + +namespace mtg { +namespace test { +namespace trie { + +using mtg::graph::traversal::Arm; +using mtg::graph::traversal::EndReason; + +inline std::string rc_of(std::string s) { ::reverse_complement(s); return s; } + +struct StringIndex { + // 2-bit packed k-mers (k <= 31 on ACGT) + using Kmer = uint64_t; + static constexpr Kmer kNone = std::numeric_limits::max(); + // packed k-mers are far from uniform in their low bits + struct KmerHash { + size_t operator()(Kmer x) const { + x ^= x >> 33; + x *= 0xff51afd7ed558ccdULL; + x ^= x >> 33; + x *= 0xc4ceb9fe1a85ec53ULL; + x ^= x >> 33; + return static_cast(x); + } + }; + using KmerSet = tsl::hopscotch_set; + + size_t k; + bool canonical; + std::map by_label; // label -> its nodes (both orientations if canonical) + + StringIndex(size_t k_, mtg::graph::DeBruijnGraph::Mode mode, + const std::vector &seqs, + const std::vector &labels) + : k(k_), canonical(mode != mtg::graph::DeBruijnGraph::BASIC) { + assert(k >= 2 && k <= 31); + for (size_t i = 0; i < seqs.size(); ++i) { + for (size_t j = 0; j + k <= seqs[i].size(); ++j) { + const std::string kmer = seqs[i].substr(j, k); + by_label[labels[i]].insert(encode(kmer)); + if (canonical) + by_label[labels[i]].insert(encode(rc_of(kmer))); + } + } + } + + static Kmer encode(const std::string &kmer) { + Kmer code = 0; + for (char c : kmer) { + const int v = c == 'A' ? 0 : c == 'C' ? 1 : c == 'G' ? 2 : c == 'T' ? 3 : -1; + if (v < 0) + return kNone; + code = (code << 2) | static_cast(v); + } + return code; + } + + bool has(const std::string &kmer) const { + const Kmer code = encode(kmer); + for (const auto &[l, kmers] : by_label) { + if (code != kNone && kmers.count(code)) + return true; + } + return false; + } + + bool carries(const std::string &label, const std::string &kmer) const { + auto it = by_label.find(label); + const Kmer code = encode(kmer); + return code != kNone && it != by_label.end() && it->second.count(code); + } + + // the labels carrying every k-mer of |seed| (the walker's derived permitted set) + std::set carriers(const std::string &seed) const { + std::set out; + for (const auto &[l, kmers] : by_label) { + bool all = seed.size() >= k; + for (size_t j = 0; all && j + k <= seed.size(); ++j) { + all = carries(l, seed.substr(j, k)); + } + if (all) + out.insert(l); + } + return out; + } + + // the seed's nodes (both orientations if canonical): re-entering one blocks a step + std::set seed_nodes(const std::string &seed) const { + std::set out; + for (size_t j = 0; j + k <= seed.size(); ++j) { + const std::string km = seed.substr(j, k); + out.insert(km); + if (canonical) + out.insert(rc_of(km)); + } + return out; + } + + std::string boundary(const std::string &seed, Arm arm) const { + return arm == Arm::RIGHT ? seed.substr(seed.size() - k) : seed.substr(0, k); + } + + // one step outward from node |u| by base |c|: the node entered and the (k+1)-mer the + // step spells + void step(Arm arm, const std::string &u, char c, std::string *v, std::string *edge) const { + if (arm == Arm::RIGHT) { + *v = u.substr(1) + c; + *edge = u + c; + } else { + *v = c + u.substr(0, k - 1); + *edge = c + u; + } + } + + // a step spelling a self-reverse-complementary (k+1)-mer; for even k also a step into + // or out of a self-reverse-complementary node + bool is_hairpin(const std::string &u, const std::string &v, const std::string &edge) const { + if (!canonical) + return false; + if (rc_of(edge) == edge) + return true; + return k % 2 == 0 && (rc_of(v) == v || rc_of(u) == u); + } + + // the edge identity (canonical (k+1)-mer in canonical regimes) and whether the step + // uses it in the flipped orientation + std::string edge_key(const std::string &edge, bool *flipped) const { + *flipped = false; + if (!canonical) + return edge; + const std::string r = rc_of(edge); + if (r < edge) { + *flipped = true; + return r; + } + return edge; + } +}; + + +inline uint8_t block_rank(EndReason r) { + switch (r) { + case EndReason::REACHED_SEED: return 3; + case EndReason::EDGE_REUSE_RC: return 2; + case EndReason::EDGE_REUSE: return 1; + default: return 0; + } +} + +// a structural successor of the current node and its admissibility under the rule +struct Step { + char ch; + std::string v; // the node entered + std::string key; // edge identity + bool flipped; + bool hairpin; // skipped + bool blocked; + EndReason block = EndReason::DEAD_END; + bool admissible() const { return !hairpin && !blocked; } +}; + +struct WalkRule { + const StringIndex &ix; + std::set seed; // the seed's nodes + Arm arm; + + WalkRule(const StringIndex &ix_, const std::string &seed_seq, Arm arm_) + : ix(ix_), seed(ix_.seed_nodes(seed_seq)), arm(arm_) {} + + // every node reachable by one step, in base order, with the rule applied against the + // edges |used| so far on this path (key -> orientation of the use) + std::vector steps(const std::string &u, const std::map &used) const { + std::vector out; + for (char c : std::string("ACGT")) { + Step s; + s.ch = c; + std::string edge; + ix.step(arm, u, c, &s.v, &edge); + if (!ix.has(s.v)) + continue; + s.hairpin = ix.is_hairpin(u, s.v, edge); + s.blocked = false; + s.key = ix.edge_key(edge, &s.flipped); + if (!s.hairpin) { + if (seed.count(s.v)) { + s.blocked = true; + s.block = EndReason::REACHED_SEED; + } else if (auto it = used.find(s.key); it != used.end()) { + s.blocked = true; + s.block = it->second != s.flipped ? EndReason::EDGE_REUSE_RC + : EndReason::EDGE_REUSE; + } + } + out.push_back(std::move(s)); + } + return out; + } +}; + +// walk (outward) -> the reason it ends +using RefTrie = std::map; + +namespace detail { + +// The end reason of a walk none of whose candidate steps is followed. |passing| are the +// steps whose node carries the label: none -> label_lost; some blocked -> the strongest +// block; only skipped hairpins -> dead_end. +inline EndReason stop_reason(const std::vector &steps, const std::vector &passing) { + bool any_pass = false; + uint8_t rank = 0; + for (size_t i = 0; i < steps.size(); ++i) { + if (!passing[i]) + continue; + any_pass = true; + if (steps[i].blocked) + rank = std::max(rank, block_rank(steps[i].block)); + } + if (!any_pass) + return EndReason::LABEL_LOST; + return rank == 3 ? EndReason::REACHED_SEED + : rank == 2 ? EndReason::EDGE_REUSE_RC + : rank == 1 ? EndReason::EDGE_REUSE : EndReason::DEAD_END; +} + +inline void dfs(const WalkRule &rule, const std::string &label, const std::string &u, + uint64_t radius, std::map &used, std::string &walk, + RefTrie *out) { + if (walk.size() >= radius) { + (*out)[walk] = EndReason::MAX_EXTENSION; + return; + } + const std::vector steps = rule.steps(u, used); + if (steps.empty()) { + (*out)[walk] = EndReason::DEAD_END; + return; + } + std::vector passing(steps.size()); + bool followed = false; + for (size_t i = 0; i < steps.size(); ++i) { + const Step &s = steps[i]; + passing[i] = rule.ix.carries(label, s.v); + if (!passing[i] || !s.admissible()) + continue; + followed = true; + used[s.key] = s.flipped; + walk.push_back(s.ch); + dfs(rule, label, s.v, radius, used, walk, out); + walk.pop_back(); + used.erase(s.key); + } + if (!followed) + (*out)[walk] = stop_reason(steps, passing); +} + +} // namespace detail + +// One label's maximal label-consistent walks, with the reason each ends. Empty when the +// label does not carry the seed boundary. +inline RefTrie label_trie(const StringIndex &ix, const std::string &label, + const std::string &seed, Arm arm, uint64_t radius) { + RefTrie out; + const std::string b = ix.boundary(seed, arm); + if (!ix.carries(label, b)) + return out; + WalkRule rule(ix, seed, arm); + std::map used; + std::string walk; + detail::dfs(rule, label, b, radius, used, walk, &out); + return out; +} + +// The claims of a permitted set: (walk, label) -> the reason the label's walk ends +using RefClaims = std::map, EndReason>; + +inline RefClaims claims(const StringIndex &ix, const std::set &permitted, + const std::string &seed, Arm arm, uint64_t radius) { + RefClaims out; + for (const std::string &l : permitted) { + for (const auto &[w, reason] : label_trie(ix, l, seed, arm, radius)) { + out[{ w, l }] = reason; + } + } + return out; +} + +} // namespace trie +} // namespace test +} // namespace mtg + +#endif // __TEST_TRIE_REFERENCE_HPP__ diff --git a/metagraph/tests/graph/traversal/test_walker.cpp b/metagraph/tests/graph/traversal/test_walker.cpp new file mode 100644 index 0000000000..9bfd59fab7 --- /dev/null +++ b/metagraph/tests/graph/traversal/test_walker.cpp @@ -0,0 +1,1049 @@ +#include "gtest/gtest.h" + +#include +#include +#include +#include +#include +#include + +#include "tests/graph/all/test_dbg_helpers.hpp" +#include "tests/annotation/test_annotated_dbg_helpers.hpp" +#include "tests/graph/traversal/walker_test_helpers.hpp" + +#include "graph/traversal/walker.hpp" +#include "graph/traversal/resolve.hpp" +#include "graph/annotated_dbg.hpp" +#include "graph/representation/canonical_dbg.hpp" +#include "graph/representation/succinct/dbg_succinct.hpp" +#include "graph/representation/hash/dbg_sshash.hpp" +#include "annotation/representation/column_compressed/annotate_column_compressed.hpp" +#include "annotation/coord_to_header.hpp" + + +namespace { + +using namespace mtg; +using namespace mtg::graph; +using namespace mtg::graph::traversal; +using namespace mtg::test; + +std::vector all_modes() { + return { + DeBruijnGraph::BASIC, +#if ! _PROTEIN_GRAPH + DeBruijnGraph::CANONICAL, + DeBruijnGraph::PRIMARY, +#endif + }; +} + +std::vector canonical_modes() { + return { +#if ! _PROTEIN_GRAPH + DeBruijnGraph::CANONICAL, + DeBruijnGraph::PRIMARY, +#endif + }; +} + +Strategy strategy(size_t max_label_branches = 0) { + Strategy st; + st.max_label_branches = max_label_branches; + return st; +} + +SeedResult run(const AnnotatedDBG &anno, + const std::string &seq, + const std::vector &labels, + const Strategy &st = Strategy()) { + LabelOracle oracle(anno); + Seed seed; + seed.sequence = seq; + seed.labels = labels; + return traverse_seed(oracle, seed, st); +} + +// every k-mer of seed + flank (or flank + seed) is in the graph +void check_spelled(const AnnotatedDBG &anno, const ArmResult &arm, const PathResult &path, + const std::string &seed) { + std::string flank = spell_path(arm, path); + EXPECT_EQ(path.length_bp, flank.size()); + std::string full = arm.arm == Arm::RIGHT ? seed + flank : flank + seed; + for (node_index n : map_to_nodes_sequentially(anno.get_graph(), full)) { + ASSERT_NE(npos, n) << full; + } +} + +size_t count_ends(const ArmResult &arm, EndReason reason) { + size_t n = 0; + for (const auto &seg : arm.segments) { + for (const auto &e : seg.label_ends) { + n += e.reason == reason; + } + } + return n; +} + +const LabelEnd* label_end(const ArmResult &arm, LabelId label) { + for (const auto &seg : arm.segments) { + for (const auto &e : seg.label_ends) { + if (e.label == label) + return &e; + } + } + return nullptr; +} + +size_t path_ends(const PathResult &path, EndReason reason) { + return std::count_if(path.end_labels.begin(), path.end_labels.end(), + [&](const PathEnd &e) { return e.reason == reason; }); +} + +// The invariants of every arm: steps are the bases added, every label end of a leaf is on +// its path, and a path's end labels are the labels at its leaf's last node +void check_invariants(const ArmResult &arm) { + uint64_t bases = 0; + for (const auto &seg : arm.segments) { + bases += seg.length_bp; + EXPECT_LE(seg.parents.size(), 1u); + for (size_t c : seg.children) { + EXPECT_EQ(std::vector{ seg.id }, arm.segments[c].parents); + } + } + EXPECT_EQ(arm.steps, bases); + EXPECT_EQ(arm.output_bp, bases); + for (const auto &path : arm.paths) { + const Segment &leaf = arm.segments[path.leaf]; + EXPECT_TRUE(leaf.children.empty()); + EXPECT_EQ(leaf.from_bp + leaf.length_bp, path.length_bp); + std::vector labels; + for (const auto &e : path.end_labels) { + labels.push_back(e.label); + const LabelEnd *end = nullptr; + for (const auto &x : leaf.label_ends) { + if (x.label == e.label && x.at_bp == path.length_bp) + end = &x; + } + ASSERT_NE(nullptr, end) << "label " << e.label; + EXPECT_EQ(end->reason, e.reason); + } + EXPECT_EQ(leaf.labels_end, labels); + } + EXPECT_EQ(arm.status == ArmResult::COMPLETE, !arm.cap_trigger.has_value()); +} + +std::string serialize(const SeedResult &r) { + std::ostringstream os; + os << r.validated_seed_id << '|' << r.label_dict.size() << '|' << r.length_bp << '|'; + for (const auto &d : r.dropped_labels) { + os << d.name << ';'; + } + for (const ArmResult &a : r.arms) { + os << "\nARM " << to_string(a.arm) << ' ' << a.requested << ' ' << a.status << ' ' + << a.steps << ' ' << a.output_bp << ' ' << a.frontier_live_paths << ' ' + << a.complete_to_bp; + if (a.cap_trigger) + os << " cap " << to_string(a.cap_trigger->reason) << ' ' << a.cap_trigger->at_bp; + for (const Segment &s : a.segments) { + os << "\n S" << s.id << " p"; + for (size_t p : s.parents) { + os << p << ','; + } + os << " c"; + for (size_t c : s.children) { + os << c << ','; + } + os << ' ' << s.from_bp << ' ' << s.length_bp << ' ' << s.sequence << " ls"; + for (LabelId l : s.labels_start) { + os << l << ','; + } + os << " le"; + for (LabelId l : s.labels_end) { + os << l << ','; + } + for (const LabelEnd &e : s.label_ends) { + os << "\n E" << e.at_bp << ' ' << e.label << ' ' << to_string(e.reason) + << ' ' << e.text; + } + } + for (const Split &sp : a.splits) { + os << "\n X" << sp.at_bp << ' ' << sp.segment << ' ' << sp.ambiguous << ' '; + for (size_t c : sp.children) { + os << c << ','; + } + } + for (const PathResult &p : a.paths) { + os << "\n P" << p.id << ' ' << p.leaf << ' ' << p.length_bp << ' ' + << (p.path_reason ? to_string(*p.path_reason) : "-") << ' '; + for (const auto &e : p.end_labels) { + os << e.label << ':' << to_string(e.reason) << ':' << e.branches << ','; + } + } + } + return os.str(); +} + + +template +class WalkerTest : public ::testing::Test {}; + +typedef ::testing::Types< + std::pair>, + std::pair>, + std::pair> +> WalkerTypes; +TYPED_TEST_SUITE(WalkerTest, WalkerTypes); + + +// one record, a seed in its middle: both arms run to the record's ends +TYPED_TEST(WalkerTest, Linear) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + std::string S = clean_block(200, 1); + std::string seed = S.substr(60, 40); + for (auto mode : all_modes()) { + auto anno = build_anno_graph(kK, { S }, { "A" }, mode); + auto res = run(*anno, seed, { "A" }); + ASSERT_EQ(1u, res.label_dict.size()); + EXPECT_TRUE(res.dropped_labels.empty()); + EXPECT_EQ(16u, res.validated_seed_id.size()); + EXPECT_FALSE(res.seed_id_mismatch); + EXPECT_EQ(40u, res.length_bp); + EXPECT_EQ(30u, res.num_kmers); + + const ArmResult &right = res.arms[kRight], &left = res.arms[kLeft]; + for (const ArmResult *arm : { &left, &right }) { + ASSERT_TRUE(arm->requested); + EXPECT_EQ(ArmResult::COMPLETE, arm->status); + EXPECT_EQ(Strategy().max_extension_bp, arm->complete_to_bp); + ASSERT_EQ(1u, arm->paths.size()); + ASSERT_EQ(1u, arm->segments.size()); + const PathResult &path = arm->paths[0]; + EXPECT_FALSE(path.path_reason.has_value()); + ASSERT_EQ(1u, path.end_labels.size()); + EXPECT_EQ(0u, path.end_labels[0].label); + EXPECT_EQ(EndReason::DEAD_END, path.end_labels[0].reason); + EXPECT_EQ(0u, path.end_labels[0].branches); + EXPECT_EQ((std::vector{ 0 }), arm->segments[0].labels_start); + EXPECT_EQ((std::vector{ 0 }), arm->segments[0].labels_end); + EXPECT_TRUE(arm->splits.empty()); + check_spelled(*anno, *arm, path, seed); + check_invariants(*arm); + EXPECT_EQ(arm->steps, path.length_bp); + } + EXPECT_EQ(S.substr(100), spell_path(right, right.paths[0])) << mode; + EXPECT_EQ(S.substr(0, 60), spell_path(left, left.paths[0])) << mode; + + // single-direction requests produce one arm only + Strategy st; + st.direction = Strategy::LEFT; + auto only_left = run(*anno, seed, { "A" }, st); + EXPECT_FALSE(only_left.arms[kRight].requested); + EXPECT_TRUE(only_left.arms[kRight].paths.empty()); + EXPECT_TRUE(only_left.arms[kRight].segments.empty()); + ASSERT_EQ(1u, only_left.arms[kLeft].paths.size()); + EXPECT_EQ(S.substr(0, 60), spell_path(only_left.arms[kLeft], only_left.arms[kLeft].paths[0])); + st.direction = Strategy::RIGHT; + auto only_right = run(*anno, seed, { "A" }, st); + EXPECT_FALSE(only_right.arms[kLeft].requested); + ASSERT_EQ(1u, only_right.arms[kRight].paths.size()); + EXPECT_EQ(S.substr(100), spell_path(only_right.arms[kRight], only_right.arms[kRight].paths[0])); + + // a lower-case seed maps the same way and gets the same seed id + std::string lower = seed; + for (char &c : lower) { + c = std::tolower(c); + } + auto lc = run(*anno, lower, { "A" }); + EXPECT_EQ(serialize(res), serialize(lc)); + EXPECT_EQ(make_seed_id(seed, mode != DeBruijnGraph::BASIC, { "A" }), lc.validated_seed_id); + LabelOracle oracle(*anno); + Seed s; + s.sequence = lower; + s.labels = { "A" }; + s.seed_id = res.validated_seed_id; + EXPECT_FALSE(traverse_seed(oracle, s, Strategy()).seed_id_mismatch); + s.seed_id = "0000000000000000"; + auto mismatch = traverse_seed(oracle, s, Strategy()); + EXPECT_TRUE(mismatch.seed_id_mismatch); + EXPECT_EQ("0000000000000000", mismatch.seed_id); + } +} + +// exact label-end positions: A on W·R·X, B on R·X·Y, seed R +TYPED_TEST(WalkerTest, LabelEndCoordinates) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + LabelEndFixture f; + for (auto mode : all_modes()) { + auto anno = build_anno_graph(kK, f.seqs(), f.labels(), mode); + auto res = run(*anno, f.R, { "A", "B" }); + const ArmResult &right = res.arms[kRight], &left = res.arms[kLeft]; + check_invariants(right); + check_invariants(left); + // right: both carry X, A ends after it, B goes on through Y + ASSERT_EQ(1u, right.paths.size()) << mode; + EXPECT_EQ(f.X + f.Y, spell_path(right, right.paths[0])); + const LabelEnd *a = label_end(right, 0), *b = label_end(right, 1); + ASSERT_TRUE(a && b); + EXPECT_EQ(f.X.size(), a->at_bp); + EXPECT_EQ(EndReason::LABEL_LOST, a->reason); + EXPECT_EQ(f.X.size() + f.Y.size(), b->at_bp); + EXPECT_EQ(EndReason::DEAD_END, b->reason); + ASSERT_EQ(1u, right.paths[0].end_labels.size()); + EXPECT_EQ(1u, right.paths[0].end_labels[0].label); + // left: only A carries W; B is lost at the seed boundary + ASSERT_EQ(1u, left.paths.size()); + EXPECT_EQ(f.W, spell_path(left, left.paths[0])); + EXPECT_EQ(0u, label_end(left, 1)->at_bp); + EXPECT_EQ(EndReason::LABEL_LOST, label_end(left, 1)->reason); + EXPECT_EQ(f.W.size(), label_end(left, 0)->at_bp); + EXPECT_EQ((std::vector{ 0, 1 }), left.segments[0].labels_start); + EXPECT_EQ((std::vector{ 0 }), left.segments[0].labels_end); + } +} + +// two labels diverging at the seed boundary: a split, one label per child +TYPED_TEST(WalkerTest, Divergence) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + auto b = fork_blocks(3); + const std::string &X = b[0], &P = b[1], &Q = b[2]; + for (auto mode : all_modes()) { + auto anno = build_anno_graph(kK, { X + P, X + Q }, { "A", "B" }, mode); + auto res = run(*anno, X, { "A", "B" }); + const ArmResult &right = res.arms[kRight]; + check_invariants(right); + ASSERT_EQ(1u, right.splits.size()); + EXPECT_EQ(0u, right.splits[0].at_bp); + EXPECT_EQ(0u, right.splits[0].segment); + EXPECT_FALSE(right.splits[0].ambiguous); + ASSERT_EQ(2u, right.splits[0].children.size()); + ASSERT_EQ(3u, right.segments.size()); + EXPECT_EQ(0u, right.segments[0].length_bp); + EXPECT_EQ((std::vector{ 0, 1 }), right.segments[0].labels_end); + ASSERT_EQ(2u, right.paths.size()); + std::map leaves; + for (const auto &path : right.paths) { + ASSERT_EQ(1u, path.end_labels.size()); + EXPECT_EQ(EndReason::DEAD_END, path.end_labels[0].reason); + EXPECT_EQ(0u, path.end_labels[0].branches); + leaves[spell_path(right, path)] = path.end_labels[0].label; + check_spelled(*anno, right, path, X); + } + ASSERT_EQ(2u, leaves.size()); + EXPECT_EQ(0u, leaves.at(P)); + EXPECT_EQ(1u, leaves.at(Q)); + EXPECT_EQ(0u, count_ends(right, EndReason::BRANCH)); + EXPECT_EQ(80u, right.steps); + } +} + +// one label on both sides of a structural branch +TYPED_TEST(WalkerTest, AmbiguousBranch) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + std::vector b; + for (uint32_t seed = 4; ; ++seed) { + b = clean_blocks({ 30, 25, 30, 25, 30 }, seed); + if (b[1][0] != b[3][0]) + break; + } + const std::string &X = b[0], &P = b[1], &Y = b[2], &Q = b[3], &Z = b[4]; + for (auto mode : all_modes()) { + auto anno = build_anno_graph(kK, { X + P + Y, X + Q + Z }, + { "A", "A" }, mode); + // limit 0: the label stops at the branch + auto res = run(*anno, X, { "A" }, strategy(0)); + const ArmResult &right = res.arms[kRight]; + check_invariants(right); + ASSERT_EQ(1u, right.paths.size()); + EXPECT_EQ(0u, right.paths[0].length_bp); + EXPECT_EQ(1u, path_ends(right.paths[0], EndReason::BRANCH)); + EXPECT_EQ(0u, label_end(right, 0)->at_bp); + EXPECT_TRUE(right.splits.empty()); + EXPECT_EQ(0u, right.steps); + + // limit 1 and unlimited: both continuations are followed, each charged one branch + for (size_t limit : { size_t(1), Strategy::kUnlimited }) { + res = run(*anno, X, { "A" }, strategy(limit)); + const ArmResult &r1 = res.arms[kRight]; + check_invariants(r1); + ASSERT_EQ(2u, r1.paths.size()); + std::set flanks; + for (const auto &path : r1.paths) { + flanks.insert(spell_path(r1, path)); + ASSERT_EQ(1u, path.end_labels.size()); + EXPECT_EQ(1u, path.end_labels[0].branches); + EXPECT_EQ(EndReason::DEAD_END, path.end_labels[0].reason); + check_spelled(*anno, r1, path, X); + } + EXPECT_EQ((std::set{ P + Y, Q + Z }), flanks); + ASSERT_EQ(1u, r1.splits.size()); + EXPECT_TRUE(r1.splits[0].ambiguous); + EXPECT_EQ(0u, count_ends(r1, EndReason::BRANCH)); + } + } +} + +// a cycle entered through junction J: the walk rule lets it be walked once +TYPED_TEST(WalkerTest, CycleJunction) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + std::vector b; + for (uint32_t seed = 6; ; ++seed) { + b = clean_blocks({ 30, kK, 40, 30 }, seed); + if (b[2][0] != b[3][0]) + break; + } + const std::string &L = b[0], &J = b[1], &C = b[2], &E = b[3]; + std::string S = L + J + C + J + E; + for (auto mode : all_modes()) { + auto anno = build_anno_graph(kK, { S }, { "A" }, mode); + auto res = run(*anno, L, { "A" }, strategy(0)); + const ArmResult &r0 = res.arms[kRight]; + check_invariants(r0); + ASSERT_EQ(1u, r0.paths.size()); + EXPECT_EQ(J, spell_path(r0, r0.paths[0])); + EXPECT_EQ(1u, count_ends(r0, EndReason::BRANCH)); + + res = run(*anno, L, { "A" }, strategy(1)); + const ArmResult &r1 = res.arms[kRight]; + check_invariants(r1); + ASSERT_EQ(2u, r1.paths.size()); + std::set flanks; + for (const auto &path : r1.paths) { + flanks.insert(spell_path(r1, path)); + ASSERT_EQ(1u, path.end_labels.size()); + EXPECT_EQ(1u, path.end_labels[0].branches); + EXPECT_EQ(EndReason::DEAD_END, path.end_labels[0].reason); + check_spelled(*anno, r1, path, L); + } + // the second time at J the cycle's first edge was taken: only the exit is followed + EXPECT_EQ((std::set{ J + E, J + C + J + E }), flanks); + EXPECT_EQ(0u, count_ends(r1, EndReason::BRANCH)); + EXPECT_EQ(0u, count_ends(r1, EndReason::EDGE_REUSE)); + } +} + +// a homopolymer self-loop: each loop step re-takes the loop edge +TYPED_TEST(WalkerTest, SelfLoop) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + std::vector b; + for (uint32_t seed = 7; ; ++seed) { + b = clean_blocks({ 30, 30 }, seed); + if (b[0].back() != 'A' && b[1][0] != 'A' && b[0].back() != 'T') + break; + } + const std::string &L = b[0], &E = b[1]; + std::string H(kK + 1, 'A'); + std::string S = L + H + E; + for (auto mode : all_modes()) { + auto anno = build_anno_graph(kK, { S }, { "A" }, mode); + auto res = run(*anno, L, { "A" }, strategy(0)); + const ArmResult &r0 = res.arms[kRight]; + check_invariants(r0); + ASSERT_EQ(1u, r0.paths.size()); + EXPECT_EQ(kK - 1, r0.paths[0].length_bp); + EXPECT_EQ(1u, count_ends(r0, EndReason::BRANCH)); + check_spelled(*anno, r0, r0.paths[0], L); + + res = run(*anno, L, { "A" }, strategy(1)); + const ArmResult &r1 = res.arms[kRight]; + check_invariants(r1); + ASSERT_EQ(2u, r1.paths.size()); + std::set flanks; + for (const auto &path : r1.paths) { + flanks.insert(spell_path(r1, path)); + check_spelled(*anno, r1, path, L); + } + EXPECT_EQ((std::set{ H.substr(1), H.substr(2) + E }), flanks); + EXPECT_EQ(1u, count_ends(r1, EndReason::BRANCH)); + EXPECT_EQ(1u, count_ends(r1, EndReason::DEAD_END)); + + res = run(*anno, L, { "A" }, strategy(2)); + const ArmResult &r2 = res.arms[kRight]; + check_invariants(r2); + ASSERT_EQ(3u, r2.paths.size()); + flanks.clear(); + for (const auto &path : r2.paths) { + flanks.insert(spell_path(r2, path)); + check_spelled(*anno, r2, path, L); + EXPECT_EQ(1u, path_ends(path, EndReason::DEAD_END)); + } + EXPECT_EQ((std::set{ H + E, H.substr(1) + E, H.substr(2) + E }), flanks); + EXPECT_EQ(0u, count_ends(r2, EndReason::BRANCH)); + } +} + +// a circular record: both arms go round and stop where they would step back into the seed +TYPED_TEST(WalkerTest, Circle) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + std::string S = clean_block(150, 8); + std::string record = S + S.substr(0, kK - 1 + 10); + std::string seed = S.substr(50, 40); + for (auto mode : all_modes()) { + auto anno = build_anno_graph(kK, { record }, { "A" }, mode); + auto res = run(*anno, seed, { "A" }); + for (size_t a : { kLeft, kRight }) { + const ArmResult &arm = res.arms[a]; + check_invariants(arm); + ASSERT_EQ(1u, arm.paths.size()) << mode << " arm " << a; + EXPECT_EQ(ArmResult::COMPLETE, arm.status); + EXPECT_EQ(1u, path_ends(arm.paths[0], EndReason::REACHED_SEED)); + EXPECT_EQ(S.size() - seed.size() + kK - 1, arm.paths[0].length_bp); + check_spelled(*anno, arm, arm.paths[0], seed); + } + EXPECT_EQ(S.substr(90) + S.substr(0, 60), spell_path(res.arms[kRight], res.arms[kRight].paths[0])); + EXPECT_EQ(S.substr(80) + S.substr(0, 50), spell_path(res.arms[kLeft], res.arms[kLeft].paths[0])); + } +} + +// a self-reverse-complementary stretch: the hairpin step back is never taken +TYPED_TEST(WalkerTest, Hairpin) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + std::string X = "AACGTACGTT"; + ASSERT_EQ(X, rc(X)); + ASSERT_EQ(kK - 1, X.size()); + std::string L, R; + for (uint32_t seed = 9; ; ++seed) { + L = clean_block(40, seed); + R = clean_block(40, seed + 100); + if (R[0] != rc(std::string(1, L.back()))[0] && is_clean(L + X.substr(0, 3)) + && is_clean(X.substr(7) + R)) + break; + } + std::string S = L + X + R; + std::string seed = L.substr(0, 30); + for (auto mode : canonical_modes()) { + auto anno = build_anno_graph(kK, { S }, { "A" }, mode); + auto res = run(*anno, seed, { "A" }); + const ArmResult &right = res.arms[kRight]; + check_invariants(right); + ASSERT_EQ(1u, right.paths.size()) << mode; + EXPECT_EQ(S.substr(30), spell_path(right, right.paths[0])) << mode; + EXPECT_EQ(0u, count_ends(right, EndReason::BRANCH)); + EXPECT_TRUE(right.splits.empty()); + } +} + +// CLI-built DBGSuccinct graphs keep their dummy k-mers: the walk is the same +TEST(Walker, UnmaskedSuccinct) { + LabelEndFixture f; + for (auto mode : all_modes()) { + auto masked = build_anno_graph>(kK, f.seqs(), f.labels(), mode); + std::shared_ptr graph; + std::shared_ptr base; + if (mode == DeBruijnGraph::PRIMARY) { + auto canonical = build_graph(kK, f.seqs(), DeBruijnGraph::CANONICAL); + std::vector contigs; + canonical->call_sequences([&](const std::string &c, const auto &) { contigs.push_back(c); }, 1, true); + base = std::make_shared(kK, DeBruijnGraph::PRIMARY); + for (const auto &c : contigs) { + base->add_sequence(c); + } + graph = std::make_shared(base); + } else { + base = std::make_shared(kK, mode); + for (const auto &s : f.seqs()) { + base->add_sequence(s); + } + graph = base; + } + ASSERT_EQ(nullptr, base->get_mask()); + auto unmasked = std::make_unique( + graph, std::make_unique>(base->max_index())); + for (size_t i = 0; i < f.seqs().size(); ++i) { + unmasked->annotate_sequence(f.seqs()[i], { f.labels()[i] }); + } + for (const std::string &seed : { f.R, f.X, f.W + f.R }) { + std::vector labels = seed == f.W + f.R ? std::vector{ "A" } + : std::vector{ "A", "B" }; + EXPECT_EQ(serialize(run(*masked, seed, labels)), serialize(run(*unmasked, seed, labels))) + << mode << " seed " << seed; + } + } +} + +TEST(Walker, InconsistentPrimaryGraph) { + std::string S = clean_block(80, 17); + auto base = std::make_shared(kK, DeBruijnGraph::PRIMARY); + base->add_sequence(S); + base->add_sequence(rc(S)); + base->mask_dummy_kmers(1, false); + auto graph = std::make_shared(base); + auto anno = std::make_unique( + graph, std::make_unique>(base->max_index())); + anno->annotate_sequence(S, { "A" }); + LabelOracle oracle(*anno); + ASSERT_EQ(Regime::PRIMARY, oracle.regime()); + Seed seed; + seed.sequence = S.substr(20, 30); + seed.labels = { "A" }; + try { + traverse_seed(oracle, seed, Strategy()); + FAIL() << "no exception"; + } catch (const std::runtime_error &e) { + std::string what = e.what(); + EXPECT_NE(std::string::npos, what.find("Inconsistent primary graph at node")) << what; + EXPECT_NE(std::string::npos, what.find("on arm")) << what; + } +} + +TYPED_TEST(WalkerTest, Caps) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + std::string S = clean_block(200, 10); + std::string seed = S.substr(60, 40); + auto b = fork_blocks(11); + const std::string &X = b[0], &P = b[1], &Q = b[2]; + for (auto mode : all_modes()) { + auto linear = build_anno_graph(kK, { S }, { "A" }, mode); + auto fork = build_anno_graph(kK, { X + P, X + Q }, { "A", "B" }, mode); + + // the extension length completes the arm + Strategy st; + st.max_extension_bp = 10; + auto res = run(*linear, seed, { "A" }, st); + for (size_t a : { kLeft, kRight }) { + const ArmResult &arm = res.arms[a]; + check_invariants(arm); + EXPECT_EQ(ArmResult::COMPLETE, arm.status); + EXPECT_EQ(10u, arm.complete_to_bp); + ASSERT_EQ(1u, arm.paths.size()); + EXPECT_EQ(10u, arm.paths[0].length_bp); + ASSERT_TRUE(arm.paths[0].path_reason.has_value()); + EXPECT_EQ(EndReason::MAX_EXTENSION, *arm.paths[0].path_reason); + EXPECT_EQ(1u, path_ends(arm.paths[0], EndReason::MAX_EXTENSION)); + EXPECT_EQ(10u, arm.steps); + } + + // the steps cap stops the seed, both arms + st = Strategy(); + st.max_steps = 5; + res = run(*linear, seed, { "A" }, st); + EXPECT_EQ(5u, res.arms[kLeft].steps + res.arms[kRight].steps); + for (size_t a : { kLeft, kRight }) { + const ArmResult &arm = res.arms[a]; + check_invariants(arm); + EXPECT_EQ(ArmResult::TRUNCATED, arm.status); + ASSERT_TRUE(arm.cap_trigger.has_value()); + EXPECT_EQ(EndReason::MAX_STEPS, arm.cap_trigger->reason); + EXPECT_EQ(arm.cap_trigger->at_bp, arm.complete_to_bp); + EXPECT_EQ(1u, arm.frontier_live_paths); + ASSERT_EQ(1u, arm.paths.size()); + EXPECT_EQ(EndReason::MAX_STEPS, *arm.paths[0].path_reason); + EXPECT_EQ(1u, count_ends(arm, EndReason::MAX_STEPS)); + } + + st = Strategy(); + st.max_output_bp = 10; + res = run(*linear, seed, { "A" }, st); + for (size_t a : { kLeft, kRight }) { + const ArmResult &arm = res.arms[a]; + check_invariants(arm); + EXPECT_EQ(ArmResult::TRUNCATED, arm.status); + EXPECT_EQ(10u, arm.output_bp); + ASSERT_TRUE(arm.cap_trigger.has_value()); + EXPECT_EQ(EndReason::MAX_OUTPUT, arm.cap_trigger->reason); + EXPECT_EQ(10u, arm.cap_trigger->at_bp); + EXPECT_EQ(1u, count_ends(arm, EndReason::MAX_OUTPUT)); + } + + // the live-path and path caps stop one arm only + st = Strategy(); + st.max_live_paths = 1; + res = run(*fork, X, { "A", "B" }, st); + { + const ArmResult &right = res.arms[kRight], &left = res.arms[kLeft]; + check_invariants(right); + EXPECT_EQ(ArmResult::TRUNCATED, right.status); + ASSERT_TRUE(right.cap_trigger.has_value()); + EXPECT_EQ(EndReason::MAX_LIVE_PATHS, right.cap_trigger->reason); + EXPECT_EQ(0u, right.cap_trigger->at_bp); + EXPECT_EQ(0u, right.complete_to_bp); + EXPECT_EQ(2u, count_ends(right, EndReason::MAX_LIVE_PATHS)); + EXPECT_TRUE(right.splits.empty()); + EXPECT_EQ(ArmResult::COMPLETE, left.status); + EXPECT_EQ(2u, count_ends(left, EndReason::DEAD_END)); + } + st = Strategy(); + st.max_paths = 1; + res = run(*fork, X, { "A", "B" }, st); + EXPECT_EQ(ArmResult::TRUNCATED, res.arms[kRight].status); + ASSERT_TRUE(res.arms[kRight].cap_trigger.has_value()); + EXPECT_EQ(EndReason::MAX_PATHS, res.arms[kRight].cap_trigger->reason); + EXPECT_EQ(2u, count_ends(res.arms[kRight], EndReason::MAX_PATHS)); + EXPECT_EQ(ArmResult::COMPLETE, res.arms[kLeft].status); + st.max_paths = 2; + res = run(*fork, X, { "A", "B" }, st); + EXPECT_EQ(ArmResult::COMPLETE, res.arms[kRight].status); + EXPECT_EQ(2u, res.arms[kRight].paths.size()); + + // a zero time budget still delivers the first level + st = Strategy(); + st.time_budget_ms = 0; + res = run(*linear, seed, { "A" }, st); + for (size_t a : { kLeft, kRight }) { + const ArmResult &arm = res.arms[a]; + EXPECT_EQ(ArmResult::TRUNCATED, arm.status); + EXPECT_EQ(1u, arm.steps); + ASSERT_TRUE(arm.cap_trigger.has_value()); + EXPECT_EQ(EndReason::TIME_BUDGET, arm.cap_trigger->reason); + EXPECT_EQ(1u, arm.cap_trigger->at_bp); + EXPECT_EQ(1u, arm.complete_to_bp); + EXPECT_EQ(1u, count_ends(arm, EndReason::TIME_BUDGET)); + } + } +} + +// a head that reached the extension length when a cap trips ends complete, not censored +TEST(Walker, HeadsAtTheRadiusAreCompleteWhenACapTrips) { + auto b = fork_blocks(62); + const std::string &X = b[0], &P = b[1], &Q = b[2]; + auto anno = build_anno_graph>( + kK, { X + P, X + Q }, { "A", "B" }, DeBruijnGraph::BASIC); + Strategy st; + st.direction = Strategy::RIGHT; + st.max_extension_bp = 2; + st.max_steps = 3; + auto res = run(*anno, X, { "A", "B" }, st); + const ArmResult &arm = res.arms[kRight]; + EXPECT_EQ(ArmResult::TRUNCATED, arm.status); + EXPECT_EQ(1u, arm.complete_to_bp); + ASSERT_EQ(2u, arm.paths.size()); + std::map ends; + for (const auto &p : arm.paths) { + ASSERT_EQ(1u, p.end_labels.size()); + ends[p.length_bp] = p.end_labels[0].reason; + } + ASSERT_EQ(2u, ends.size()); + EXPECT_STREQ("max_extension_bp", to_string(ends.at(2))); + EXPECT_STREQ("max_steps", to_string(ends.at(1))); + EXPECT_EQ(1u, arm.frontier_live_paths); + ASSERT_TRUE(arm.cap_trigger.has_value()); + EXPECT_EQ(1u, arm.cap_trigger->at_bp); +} + +// the split allowance of a path: past it the labels end with BRANCH, "split_limit" +TYPED_TEST(WalkerTest, SplitLimit) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + auto b = fork_blocks(3); + const std::string &X = b[0], &P = b[1], &Q = b[2]; + for (auto mode : all_modes()) { + auto anno = build_anno_graph(kK, { X + P, X + Q }, { "A", "B" }, mode); + Strategy st; + st.max_splits_per_path = 0; + auto res = run(*anno, X, { "A", "B" }, st); + const ArmResult &right = res.arms[kRight]; + check_invariants(right); + ASSERT_EQ(1u, right.paths.size()) << mode; + EXPECT_EQ(0u, right.paths[0].length_bp); + EXPECT_EQ(2u, path_ends(right.paths[0], EndReason::BRANCH)); + EXPECT_TRUE(right.splits.empty()); + EXPECT_EQ(0u, right.steps); + for (LabelId l : { 0u, 1u }) { + const LabelEnd *e = label_end(right, l); + ASSERT_NE(nullptr, e); + EXPECT_EQ(0u, e->at_bp); + EXPECT_EQ(EndReason::BRANCH, e->reason); + EXPECT_STREQ("split_limit", e->text); + } + st.max_splits_per_path = 1; + EXPECT_EQ(2u, run(*anno, X, { "A", "B" }, st).arms[kRight].paths.size()); + } +} + +TYPED_TEST(WalkerTest, SeedValidation) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + LabelEndFixture f; + for (auto mode : all_modes()) { + auto anno = build_anno_graph(kK, f.seqs(), f.labels(), mode); + // A does not carry the seed's last k-mers: dropped, with the k-mers it carries + std::string seed = f.R + f.X + f.Y.substr(0, 5); + auto res = run(*anno, seed, { "A", "B" }); + ASSERT_EQ(1u, res.dropped_labels.size()); + EXPECT_EQ("A", res.dropped_labels[0].name); + uint64_t supported = f.R.size() + f.X.size() - kK + 1; + EXPECT_EQ((std::vector>{ { 0, supported } }), + res.dropped_labels[0].runs); + ASSERT_EQ(1u, res.label_dict.size()); + EXPECT_EQ("B", res.label_dict[0].name); + EXPECT_EQ(1u, res.labels_supporting_total); + EXPECT_EQ(make_seed_id(seed, mode != DeBruijnGraph::BASIC, { "B" }), res.validated_seed_id); + EXPECT_EQ(f.Y.substr(5), spell_path(res.arms[kRight], res.arms[kRight].paths[0])); + + EXPECT_THROW(run(*anno, f.W + f.R, { "B" }), std::invalid_argument); // unsupported + EXPECT_THROW(run(*anno, random_seq(40, 99), { "A" }), SeedNotInGraph); + EXPECT_THROW(run(*anno, f.R + "ACGTACGTACGTACGT", { "A", "B" }), SeedNotInGraph); + EXPECT_THROW(run(*anno, f.R, { "C" }), std::invalid_argument); // unknown + EXPECT_THROW(run(*anno, f.R, { "A", "C" }), std::invalid_argument); + EXPECT_THROW(run(*anno, f.R.substr(0, kK - 1), { "A" }), std::invalid_argument); + EXPECT_THROW(run(*anno, f.R.substr(0, 20) + "N" + f.R.substr(21, 20), { "A" }), + std::invalid_argument); + EXPECT_THROW(run(*anno, "", { "A" }), std::invalid_argument); + EXPECT_THROW(run(*anno, f.R, { "A", "A" }), std::invalid_argument); // duplicate + EXPECT_EQ(2u, run(*anno, f.R, {}).label_dict.size()); + EXPECT_THROW(run(*anno, f.W + f.R + f.X + f.Y, {}), SeedDerivationError); + } +} + +// labels derived from the seed: its carriers, resubmittable as an explicit list +TYPED_TEST(WalkerTest, DeriveSeedLabels) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + DerivedFixture f; + for (auto mode : all_modes()) { + auto anno = build_anno_graph(kK, f.seqs(), f.labels(), mode); + Strategy st; + auto derived = run(*anno, f.R, {}, st); + EXPECT_TRUE(derived.labels_from_seed); + EXPECT_EQ(2u, derived.label_dict.size()) << mode; + EXPECT_EQ(2u, derived.labels_supporting_total); + EXPECT_EQ(0u, derived.labels_dropped); + EXPECT_TRUE(derived.labels_dropped_digest.empty()); + EXPECT_TRUE(derived.dropped_labels.empty()); + + std::vector names; + for (const auto &l : derived.label_dict) { + EXPECT_EQ(LabelKind::COLUMN, l.kind); // no CoordToHeader: column labels + names.push_back(l.name); + } + EXPECT_EQ((std::set{ "A", "B" }), std::set(names.begin(), names.end())); + for (size_t i = 1; i < derived.label_dict.size(); ++i) { // ascending column id + EXPECT_LT(derived.label_dict[i - 1].column, derived.label_dict[i].column); + } + + auto named = run(*anno, f.R, names, st); + EXPECT_FALSE(named.labels_from_seed); + EXPECT_EQ(serialize(named), serialize(derived)) << mode; + + auto all_three = run(*anno, f.prefix(), {}, st); + EXPECT_EQ(3u, all_three.label_dict.size()) << mode; + EXPECT_EQ(3u, all_three.labels_supporting_total); + + st.max_seed_labels = 2; + auto capped = run(*anno, f.prefix(), {}, st); + ASSERT_EQ(2u, capped.label_dict.size()) << mode; + EXPECT_EQ(3u, capped.labels_supporting_total); + EXPECT_EQ(1u, capped.labels_dropped); + EXPECT_TRUE(capped.dropped_labels.empty()); // capped, not unsupported + EXPECT_EQ(capped.label_dict[0].name, all_three.label_dict[0].name); + EXPECT_EQ(capped.label_dict[1].name, all_three.label_dict[1].name); + const std::string &cut = all_three.label_dict[2].name; + EXPECT_EQ(hex64(fnv1a64(cut + "\n")), capped.labels_dropped_digest) << mode; + st.max_seed_labels = 0; // a cap of zero would permit nothing + EXPECT_THROW(run(*anno, f.prefix(), {}, st), std::invalid_argument); + + try { + run(*anno, f.W + f.R + f.X + f.Y, {}); + FAIL() << "no carrier"; + } catch (const SeedDerivationError &e) { + EXPECT_EQ(SeedDerivationError::NO_CARRIER, e.cause()); + } + } +} + +TYPED_TEST(WalkerTest, Determinism) { + using Graph = typename TypeParam::first_type; + using Annotation = typename TypeParam::second_type; + auto b = fork_blocks(5); + const std::string &X = b[0], &P = b[1], &Q = b[2]; + for (auto mode : all_modes()) { + auto anno = build_anno_graph(kK, { X + P, X + Q, X + P + Q }, + { "A", "B", "C" }, mode); + const std::string first = serialize(run(*anno, X, { "A", "B", "C" }, strategy(1))); + for (int i = 0; i < 2; ++i) { + EXPECT_EQ(first, serialize(run(*anno, X, { "A", "B", "C" }, strategy(1)))); + } + // the label order of the request changes the ids, not the walk + auto a = run(*anno, X, { "A", "B", "C" }, strategy(1)); + auto c = run(*anno, X, { "C", "B", "A" }, strategy(1)); + EXPECT_EQ(a.validated_seed_id, c.validated_seed_id); + EXPECT_EQ(a.arms[kRight].paths.size(), c.arms[kRight].paths.size()); + EXPECT_EQ(a.arms[kRight].steps, c.arms[kRight].steps); + } +} + +TEST(WalkerDerive, ManyDerivedLabelsMatchTheExplicitList) { + constexpr size_t kN = 24; + std::vector lengths { 30 }; + lengths.insert(lengths.end(), kN, 20); + auto b = clean_blocks(lengths, 97); + const std::string S = b[0]; + + std::vector seqs, labels; + for (size_t i = 0; i < kN; ++i) { + seqs.push_back(S + b[1 + i]); + labels.push_back("full" + std::to_string(i)); + } + seqs.push_back(S.substr(0, S.size() - 10)); + labels.push_back("partial"); + auto anno = build_anno_graph>(kK, seqs, labels); + + auto derived = run(*anno, S, {}, strategy(Strategy::kUnlimited)); + EXPECT_TRUE(derived.labels_from_seed); + ASSERT_EQ(kN, derived.label_dict.size()); + EXPECT_EQ(kN, derived.labels_supporting_total); + EXPECT_EQ(kN, derived.arms[kRight].paths.size()); + + std::vector names; + for (const auto &l : derived.label_dict) { + names.push_back(l.name); + } + auto named = run(*anno, S, names, strategy(Strategy::kUnlimited)); + EXPECT_EQ(serialize(named), serialize(derived)); + + names.push_back("partial"); + auto mixed = run(*anno, S, names, strategy(Strategy::kUnlimited)); + EXPECT_EQ(kN, mixed.label_dict.size()); + ASSERT_EQ(1u, mixed.dropped_labels.size()); + EXPECT_EQ("partial", mixed.dropped_labels[0].name); + EXPECT_FALSE(mixed.dropped_labels[0].runs.empty()); +} + +TEST(WalkerDerive, DerivationHonoursTheTimeBudget) { + DerivedFixture f; + auto anno = build_anno_graph>( + kK, f.seqs(), f.labels(), DeBruijnGraph::BASIC); + Strategy st; + st.time_budget_ms = 1e-9; // spent by the time the first seed k-mer is consumed + try { + run(*anno, f.R, {}, st); + FAIL() << "the derivation ignored its time budget"; + } catch (const SeedDerivationError &e) { + EXPECT_EQ(SeedDerivationError::TIME_BUDGET, e.cause()); + EXPECT_EQ(1.0, e.observed()); + } + // explicit labels need no derivation: the walk stops at its first level instead + SeedResult named = run(*anno, f.R, { "A", "B" }, st); + EXPECT_EQ(2u, named.label_dict.size()); + EXPECT_EQ(ArmResult::TRUNCATED, named.arms[kRight].status); +} + +// a seed starting in a repetitive k-mer: the intersection starts from the narrowest row +TEST(WalkerDerive, IntersectionStartsFromTheCheapestRowOfTheFirstWindow) { + constexpr size_t kDecoys = 100; + auto b = clean_blocks({ kK, 30 }, 137); + const std::string U = b[0], V = b[1]; // U is exactly one k-mer + const std::string seed_seq = U + V; + + std::vector seqs, labels; + for (size_t i = 0; i < kDecoys; ++i) { + seqs.push_back(U); + labels.push_back("decoy" + std::to_string(i)); + } + seqs.push_back(seed_seq); + labels.push_back("carrierA"); + seqs.push_back(seed_seq); + labels.push_back("carrierB"); + auto anno = build_anno_graph>( + kK, seqs, labels, DeBruijnGraph::BASIC, true); + + const auto &encoder = anno->get_annotator().get_label_encoder(); + std::vector> headers(labels.size()); + std::vector> num_kmers(labels.size()); + for (size_t i = 0; i < labels.size(); ++i) { + Column c = static_cast(encoder.encode(labels[i])); + headers[c] = { "h_" + labels[i] }; + num_kmers[c] = { seqs[i].size() - kK + 1 }; + } + annot::CoordToHeader cth(std::move(headers), std::move(num_kmers)); + + LabelOracle oracle(*anno, &cth); + Seed seed; + seed.sequence = seed_seq; + auto derived = traverse_seed(oracle, seed, Strategy()); + ASSERT_EQ(2u, derived.label_dict.size()); + std::set names; + for (const auto &l : derived.label_dict) { + EXPECT_EQ(LabelKind::HEADER, l.kind); + names.insert(l.name); + } + EXPECT_EQ((std::set{ "h_carrierA", "h_carrierB" }), names); + // the decoys' coordinates of U were never mapped to their sequences + EXPECT_LT(derived.annotation_counters.coords_mapped, kDecoys); +} + +// a header in two columns cannot be resubmitted by name: refused, with the lever +TEST(WalkerDerive, AmbiguousDerivedHeaderNamesAreRefused) { + auto b = clean_blocks({ 30, 20 }, 211); + const std::string S = b[0], T = b[1]; + auto anno = build_anno_graph>( + kK, { S + T, S }, { "colA", "colB" }, DeBruijnGraph::BASIC, true); + + const auto &encoder = anno->get_annotator().get_label_encoder(); + std::vector> headers(2); + std::vector> num_kmers(2); + headers[static_cast(encoder.encode("colA"))] = { "ACC1" }; + num_kmers[static_cast(encoder.encode("colA"))] = { (S + T).size() - kK + 1 }; + headers[static_cast(encoder.encode("colB"))] = { "ACC1" }; + num_kmers[static_cast(encoder.encode("colB"))] = { S.size() - kK + 1 }; + annot::CoordToHeader cth(std::move(headers), std::move(num_kmers)); + + LabelOracle oracle(*anno, &cth); + Seed seed; + seed.sequence = S; // carried in full by both columns + Strategy st; + try { + traverse_seed(oracle, seed, st); + FAIL() << "expected an ambiguous derived header set to be refused"; + } catch (const SeedDerivationError &e) { + EXPECT_NE(nullptr, std::strstr(e.what(), "ACC1")) << e.what(); + EXPECT_NE(nullptr, std::strstr(e.what(), "seed_label_kind")) << e.what(); + EXPECT_EQ(SeedDerivationError::AMBIGUOUS_HEADER, e.cause()); + EXPECT_EQ("ACC1", e.subject()); + } + st.seed_label_kind = LabelKind::COLUMN; + auto by_column = traverse_seed(oracle, seed, st); + ASSERT_EQ(2u, by_column.label_dict.size()); + EXPECT_EQ("colA", by_column.label_dict[0].name); + EXPECT_EQ("colB", by_column.label_dict[1].name); + seed.labels = { "ACC1" }; + auto one = traverse_seed(oracle, seed, Strategy()); + EXPECT_EQ(1u, one.label_dict.size()); +} + +// a derived header that an explicit list would resolve to a column is refused +TEST(WalkerDerive, DerivedHeaderMustNotBeAColumnName) { + auto b = clean_blocks({ 40, 40, 30 }, 941); + const std::string &D = b[0], &M = b[1], &T = b[2]; + const std::string in_acc = D, in_b = M + T; // the seed M is only in column B's record + auto anno = build_anno_graph>( + kK, { in_acc, in_b }, { "ACC", "B" }, DeBruijnGraph::BASIC, true, { 0, 0 }); + const auto &enc = anno->get_annotator().get_label_encoder(); + std::vector> headers(2); + std::vector> num_kmers(2); + headers[enc.encode("ACC")] = { "decoy" }; + num_kmers[enc.encode("ACC")] = { in_acc.size() - kK + 1 }; + headers[enc.encode("B")] = { "ACC" }; + num_kmers[enc.encode("B")] = { in_b.size() - kK + 1 }; + annot::CoordToHeader cth(std::move(headers), std::move(num_kmers)); + LabelOracle oracle(*anno, &cth); + + Seed seed; + seed.sequence = M; + Strategy st; + st.direction = Strategy::RIGHT; + try { + traverse_seed(oracle, seed, st); + FAIL() << "a derived header that an explicit list resolves to a column was accepted"; + } catch (const SeedDerivationError &e) { + EXPECT_NE(std::string::npos, std::string(e.what()).find("'ACC'")) << e.what(); + EXPECT_NE(std::string::npos, std::string(e.what()).find("name of an annotation column")) << e.what(); + EXPECT_EQ(SeedDerivationError::AMBIGUOUS_HEADER, e.cause()); + } + st.seed_label_kind = LabelKind::COLUMN; + auto by_column = traverse_seed(oracle, seed, st); + ASSERT_EQ(1u, by_column.label_dict.size()); + EXPECT_EQ("B", by_column.label_dict[0].name); + EXPECT_EQ(T.size(), by_column.arms[kRight].paths[0].length_bp); +} + +} // namespace diff --git a/metagraph/tests/graph/traversal/test_walker_reference.cpp b/metagraph/tests/graph/traversal/test_walker_reference.cpp new file mode 100644 index 0000000000..e4f84ad65d --- /dev/null +++ b/metagraph/tests/graph/traversal/test_walker_reference.cpp @@ -0,0 +1,339 @@ +#include "gtest/gtest.h" + +#include +#include + +#include "tests/graph/all/test_dbg_helpers.hpp" +#include "tests/annotation/test_annotated_dbg_helpers.hpp" +#include "tests/graph/traversal/walker_test_helpers.hpp" +#include "tests/graph/traversal/test_trie_reference.hpp" + +#include "graph/traversal/walker.hpp" +#include "graph/annotated_dbg.hpp" +#include "annotation/representation/column_compressed/annotate_column_compressed.hpp" +#include "annotation/representation/annotation_matrix/static_annotators_def.hpp" + +/** + * The walker against an independent model (test_trie_reference.hpp): for every fixture, + * mode and permitted set, the labeled walks the walker claims with an unlimited branch + * allowance — every (walk, label) where a label ends, with its reason — are exactly the + * maximal label-consistent walks the records admit under the walk rule. And on top: + * - a run over several labels is the union of the single-label runs; + * - a permitted set derived from the seed is its carriers and runs like the explicit list; + * - with a branch allowance of 0 every claim is a claim of the reference cut at the + * label's first fork, where it ends with BRANCH. + */ + +namespace { + +using namespace mtg; +using namespace mtg::graph; +using namespace mtg::graph::traversal; +using namespace mtg::test; + +std::vector all_modes() { + return { + DeBruijnGraph::BASIC, +#if ! _PROTEIN_GRAPH + DeBruijnGraph::CANONICAL, + DeBruijnGraph::PRIMARY, +#endif + }; +} + +std::string reversed(std::string s) { + std::reverse(s.begin(), s.end()); + return s; +} + +struct Case { + std::string name; + size_t k; + std::vector sequences; + std::vector labels; + std::string seed; + uint64_t radius = 200; +}; + +// the walk (outward from the seed boundary) of |segment|'s chain up to |at_bp| +std::string outward_walk(const ArmResult &arm, size_t segment, uint64_t at_bp) { + std::vector chain; + for (size_t s = segment; ; s = arm.segments[s].parents[0]) { + chain.push_back(s); + if (arm.segments[s].parents.empty()) + break; + } + std::string walk; + for (auto it = chain.rbegin(); it != chain.rend(); ++it) { + const std::string &seq = arm.segments[*it].sequence; + walk += arm.arm == Arm::RIGHT ? seq : reversed(seq); + } + return walk.substr(0, at_bp); +} + +trie::RefClaims walker_claims(const SeedResult &res, Arm arm) { + trie::RefClaims out; + const ArmResult &a = res.arms[static_cast(arm)]; + for (const Segment &seg : a.segments) { + for (const LabelEnd &e : seg.label_ends) { + auto [it, inserted] = out.emplace( + std::make_pair(outward_walk(a, seg.id, e.at_bp), res.label_dict[e.label].name), + e.reason); + EXPECT_TRUE(inserted) << "a (walk, label) claimed twice"; + } + } + return out; +} + +std::string show(const trie::RefClaims &claims) { + std::string out; + for (const auto &[key, reason] : claims) { + out += "\n " + key.second + " " + to_string(reason) + " [" + key.first + "]"; + } + return out; +} + +SeedResult walk(const AnnotatedDBG &anno, const Case &c, const std::vector &labels, + size_t branches = Strategy::kUnlimited) { + LabelOracle oracle(anno); + Seed seed; + seed.sequence = c.seed; + seed.labels = labels; + Strategy st; + st.max_label_branches = branches; + st.max_extension_bp = c.radius; + // no cap may trip: the comparison is with the whole reference trie + st.max_splits_per_path = std::numeric_limits::max(); + st.max_live_paths = 1'000'000; + st.max_paths = 1'000'000; + return traverse_seed(oracle, seed, st); +} + +template +void check_case(const Case &c, std::vector modes = all_modes()) { + for (auto mode : modes) { + const std::string where = c.name + " mode " + std::to_string(mode); + auto anno = build_anno_graph(c.k, c.sequences, c.labels, mode); + trie::StringIndex ix(c.k, mode, c.sequences, c.labels); + const std::set carriers = ix.carriers(c.seed); + ASSERT_FALSE(carriers.empty()) << where; + const std::vector permitted(carriers.begin(), carriers.end()); + + const SeedResult all = walk(*anno, c, permitted); + for (Arm arm : { Arm::LEFT, Arm::RIGHT }) { + const trie::RefClaims expected = trie::claims(ix, carriers, c.seed, arm, c.radius); + const trie::RefClaims got = walker_claims(all, arm); + EXPECT_EQ(expected, got) << where << " arm " << to_string(arm) + << "\nexpected:" << show(expected) << "\ngot:" << show(got); + EXPECT_EQ(ArmResult::COMPLETE, all.arms[static_cast(arm)].status) << where; + + // the union of the single-label runs + trie::RefClaims single; + for (const std::string &l : permitted) { + for (const auto &claim : walker_claims(walk(*anno, c, { l }), arm)) { + single.insert(claim); + } + } + EXPECT_EQ(got, single) << where << " arm " << to_string(arm); + + // limit 0: a label's claim ends at its first fork (BRANCH) or as the reference + for (const auto &[key, reason] : walker_claims(walk(*anno, c, permitted, 0), arm)) { + const auto &[w, l] = key; + if (reason == EndReason::BRANCH) { + size_t continuing = 0; + for (const auto &[ref_key, ref_reason] : expected) { + continuing += ref_key.second == l && ref_key.first.size() > w.size() + && ref_key.first.compare(0, w.size(), w) == 0; + } + EXPECT_GE(continuing, 2u) << where << ": " << l << " [" << w << "]"; + } else { + auto it = expected.find(key); + ASSERT_NE(expected.end(), it) << where << ": " << l << " [" << w << "]"; + EXPECT_EQ(it->second, reason) << where; + } + } + } + + // a derived set is the carriers and runs as the explicit list + const SeedResult derived = walk(*anno, c, {}); + std::set names; + for (const LabelRef &l : derived.label_dict) { + names.insert(l.name); + } + EXPECT_EQ(carriers, names) << where; + for (Arm arm : { Arm::LEFT, Arm::RIGHT }) { + EXPECT_EQ(walker_claims(all, arm), walker_claims(derived, arm)) << where; + } + } +} + +template +class WalkerReference : public ::testing::Test {}; + +typedef ::testing::Types< + std::pair>, + std::pair +> WalkerReferenceTypes; +TYPED_TEST_SUITE(WalkerReference, WalkerReferenceTypes); + +#define CASE_TYPES \ + using Graph = typename TypeParam::first_type; \ + using Annotation = typename TypeParam::second_type + +template +std::vector blocks_where(const std::vector &lengths, uint32_t seed, + const Cond &cond) { + for (uint32_t s = seed; ; ++s) { + auto b = clean_blocks(lengths, s); + if (cond(b)) + return b; + } +} + +TYPED_TEST(WalkerReference, Linear) { + CASE_TYPES; + auto b = clean_blocks({ 60, 30, 80 }, 101); + check_case({ "Linear", kK, { b[0] + b[1] + b[2] }, { "A" }, b[1] }); +} + +TYPED_TEST(WalkerReference, SeedAtRecordEnds) { + CASE_TYPES; + auto b = clean_blocks({ 50, 30, 50 }, 102); + const std::string &L = b[0], &S = b[1], &R = b[2]; + check_case({ "SeedAtRecordStart", kK, { S + R }, { "A" }, S }); + check_case({ "SeedAtRecordEnd", kK, { L + S }, { "A" }, S }); + check_case({ "SeedIsRecord", kK, { S }, { "A" }, S }); +} + +TYPED_TEST(WalkerReference, RadiusAtEveryBoundary) { + CASE_TYPES; + auto b = clean_blocks({ 40, 30, 50 }, 103); + const std::string &L = b[0], &S = b[1], &R = b[2]; + for (uint64_t radius : { 0ul, 1ul, 2ul, R.size() - 1, R.size(), R.size() + 1 }) { + check_case({ "Radius" + std::to_string(radius), kK, + { L + S + R }, { "A" }, S, radius }); + } +} + +TYPED_TEST(WalkerReference, Fork) { + CASE_TYPES; + auto b = fork_blocks(104); + const std::string &X = b[0], &P = b[1], &Q = b[2]; + check_case({ "Fork", kK, { X + P, X + Q }, { "A", "B" }, X }); + check_case({ "ForkOneLabel", kK, { X + P, X + Q }, { "A", "A" }, X }); +} + +TYPED_TEST(WalkerReference, Bubbles) { + CASE_TYPES; + auto b = blocks_where({ 30, 15, 22, 30, 30 }, 105, [](const auto &b) { return b[1][0] != b[2][0]; }); + const std::string &X = b[0], &P = b[1], &Q = b[2], &Y = b[3], &Z = b[4]; + check_case({ "UnequalBubble", kK, { X + P + Y + Z, X + Q + Y + Z }, + { "A", "B" }, X }); + check_case({ "UnequalBubbleOneLabel", kK, { X + P + Y + Z, X + Q + Y + Z }, + { "A", "A" }, X }); + auto n = blocks_where({ 30, 20, 20, 20, 25, 30 }, 106, [](const auto &b) { + return b[1][0] != b[4][0] && b[2][0] != b[3][0]; + }); + check_case({ "NestedBubbles", kK, + { n[0] + n[1] + n[2] + n[5], n[0] + n[1] + n[3] + n[5], + n[0] + n[4] + n[5] }, + { "A", "B", "C" }, n[0] }); +} + +TYPED_TEST(WalkerReference, TipAtTheSeedBoundary) { + CASE_TYPES; + auto b = clean_blocks({ 40, 40 }, 107); + const std::string &X = b[0], &R = b[1]; + char t = 'A'; + while (t == R[0]) { + t = "ACGT"[(std::string("ACGT").find(t) + 1) % 4]; + } + check_case({ "Tip", kK, { X + R, X + t }, { "A", "B" }, X }); +} + +TYPED_TEST(WalkerReference, HomopolymerSelfLoop) { + CASE_TYPES; + auto b = blocks_where({ 30, 30 }, 108, [](const auto &b) { + return b[0].back() != 'A' && b[1][0] != 'A' && b[0].back() != 'T' && b[1][0] != 'T'; + }); + check_case({ "HomopolymerSelfLoop", kK, + { b[0] + std::string(kK + 3, 'A') + b[1] }, { "A" }, b[0] }); +} + +TYPED_TEST(WalkerReference, CycleJunction) { + CASE_TYPES; + auto b = blocks_where({ 30, kK, 40, 30 }, 109, [](const auto &b) { return b[2][0] != b[3][0]; }); + const std::string &L = b[0], &J = b[1], &C = b[2], &E = b[3]; + check_case({ "CycleJunction", kK, { L + J + C + J + E }, { "A" }, L }); +} + +TYPED_TEST(WalkerReference, Circle) { + CASE_TYPES; + std::string S = clean_block(150, 110); + check_case({ "Circle", kK, { S + S.substr(0, kK - 1 + 10) }, { "A" }, + S.substr(50, 40) }); +} + +TYPED_TEST(WalkerReference, Repeats) { + CASE_TYPES; + auto b = blocks_where({ 30, 30, 30, 30, 30 }, 111, [](const auto &b) { + return b[2][0] != b[4][0] && b[0].back() != b[3].back(); + }); + const std::string &X = b[0], &R = b[1], &Y = b[2], &Z = b[3], &W = b[4]; + check_case({ "RepeatInTwoRecordsAndAChimera", kK, + { X + R + Y, Z + R + W, X + R + W }, { "A", "B", "C" }, R }); + auto f = fork_blocks(112); + check_case({ "RepeatedSeedInOneRecord", kK, + { f[0] + f[1] + f[0] + f[2] }, { "A" }, f[0] }); +} + +TYPED_TEST(WalkerReference, ReverseComplementRecord) { + CASE_TYPES; + auto b = fork_blocks(115); + const std::string &X = b[0], &P = b[1], &Q = b[2]; + check_case({ "ReverseComplementRecord", kK, { X + P, rc(X + Q) }, + { "A", "B" }, X }); +} + +TYPED_TEST(WalkerReference, Hairpin) { + CASE_TYPES; + const std::string X = "AACGTACGTT"; + ASSERT_EQ(X, rc(X)); + std::string L, R; + for (uint32_t seed = 116; ; ++seed) { + L = clean_block(40, seed); + R = clean_block(40, seed + 100); + if (R[0] != rc(std::string(1, L.back()))[0] && is_clean(L + X.substr(0, 3)) + && is_clean(X.substr(7) + R)) + break; + } + check_case({ "Hairpin", kK, { L + X + R }, { "A" }, L.substr(0, 30) }); +} + +// random records at k = 7 over a small alphabet of shared pieces: repeats, cycles, both +// strands and forks everywhere +TYPED_TEST(WalkerReference, RandomGraphs) { + CASE_TYPES; + std::mt19937 gen(118); + for (size_t trial = 0; trial < 12; ++trial) { + std::vector pieces; + for (size_t i = 0; i < 6; ++i) { + pieces.push_back(random_seq(6 + gen() % 10, gen())); + } + std::vector seqs, labels; + for (size_t r = 0; r < 4; ++r) { + std::string s; + for (size_t j = 0; j < 5; ++j) { + s += pieces[gen() % pieces.size()]; + } + seqs.push_back(s); + labels.push_back(std::string(1, 'A' + gen() % 3)); + } + const std::string &host = seqs[gen() % seqs.size()]; + const size_t at = gen() % (host.size() - 10); + check_case({ "Random" + std::to_string(trial), 7, seqs, labels, + host.substr(at, 10), 25 }); + } +} + +} // namespace diff --git a/metagraph/tests/graph/traversal/walker_test_helpers.hpp b/metagraph/tests/graph/traversal/walker_test_helpers.hpp new file mode 100644 index 0000000000..82d9bdbb2b --- /dev/null +++ b/metagraph/tests/graph/traversal/walker_test_helpers.hpp @@ -0,0 +1,132 @@ +#ifndef __TESTS_GRAPH_TRAVERSAL_WALKER_TEST_HELPERS_HPP__ +#define __TESTS_GRAPH_TRAVERSAL_WALKER_TEST_HELPERS_HPP__ + +#include +#include +#include +#include +#include + +#include "graph/traversal/walker.hpp" +#include "common/seq_tools/reverse_complement.hpp" + + +namespace mtg { +namespace test { + +constexpr size_t kK = 11; +constexpr size_t kLeft = static_cast(graph::traversal::Arm::LEFT); +constexpr size_t kRight = static_cast(graph::traversal::Arm::RIGHT); + +inline std::string rc(std::string s) { + ::reverse_complement(s); + return s; +} + +inline std::string random_seq(size_t len, uint32_t seed) { + std::mt19937 gen(seed); + std::string s(len, 'A'); + for (char &c : s) { + c = "ACGT"[gen() % 4]; + } + return s; +} + +// No repeated (k-1)-mer or k-mer in either orientation and no reverse-complement +// palindromes of length k-1 or k+1: every k-mer is one node with one successor in every +// graph mode, and no hairpin +inline bool is_clean(const std::string &s) { + for (size_t len : { kK - 1, kK }) { + std::set seen; + for (size_t i = 0; i + len <= s.size(); ++i) { + std::string km = s.substr(i, len); + if (seen.count(km) || seen.count(rc(km))) + return false; + seen.insert(km); + } + } + for (size_t len : { kK - 1, kK + 1 }) { + for (size_t i = 0; i + len <= s.size(); ++i) { + std::string x = s.substr(i, len); + if (rc(x) == x) + return false; + } + } + return true; +} + +inline std::string clean_block(size_t len, uint32_t seed) { + for (uint32_t s = seed; ; ++s) { + std::string b = random_seq(len, s); + if (is_clean(b)) + return b; + } +} + +// blocks of the given lengths cut from one clean sequence: clean together, in any order +// that keeps their junctions apart +inline std::vector clean_blocks(const std::vector &lengths, uint32_t seed) { + size_t total = 0; + for (size_t l : lengths) { + total += l; + } + std::string master = clean_block(total, seed); + std::vector blocks; + size_t pos = 0; + for (size_t l : lengths) { + blocks.push_back(master.substr(pos, l)); + pos += l; + } + return blocks; +} + +// X, P, Q with P and Q starting with different bases: X + P and X + Q fork after X +inline std::vector fork_blocks(uint32_t seed, + std::vector lengths = { 30, 40, 40 }) { + for (uint32_t s = seed; ; ++s) { + auto b = clean_blocks(lengths, s); + if (b[1][0] != b[2][0]) + return b; + } +} + +// A on W·R·X, B on R·X·Y +struct LabelEndFixture { + std::string W, R, X, Y; + LabelEndFixture() { + auto b = clean_blocks({ 20, 60, 25, 30 }, 2); + W = b[0]; R = b[1]; X = b[2]; Y = b[3]; + } + std::vector seqs() const { return { W + R + X, R + X + Y }; } + std::vector labels() const { return { "A", "B" }; } +}; + +// ... and C on the first 40 bases of R +struct DerivedFixture : public LabelEndFixture { + std::string prefix() const { return R.substr(0, 40); } + std::vector seqs() const { return { W + R + X, R + X + Y, prefix() }; } + std::vector labels() const { return { "A", "B", "C" }; } +}; + +// the bases of a path from the seed boundary outward (natural orientation on the left arm) +inline std::string spell_path(const graph::traversal::ArmResult &arm, + const graph::traversal::PathResult &path) { + std::vector pieces; + for (size_t s = path.leaf; ; s = arm.segments[s].parents[0]) { + pieces.push_back(&arm.segments[s].sequence); + if (arm.segments[s].parents.empty()) + break; + } + if (arm.arm == graph::traversal::Arm::RIGHT) + std::reverse(pieces.begin(), pieces.end()); + std::string out; + for (const std::string *p : pieces) { + out += *p; + } + return out; +} + +} // namespace test +} // namespace mtg + +#endif // __TESTS_GRAPH_TRAVERSAL_WALKER_TEST_HELPERS_HPP__ From cf6790ead45ee821553f034832961432666f4be9 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sun, 1 Nov 2026 23:25:35 +0100 Subject: [PATCH 03/52] Rebuild the label encoder's hash buckets on load MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit An LE-v2.0 label encoder is deserialized with the hash buckets its writer stored, which are only valid for the writer's std::hash. An annotation written by another toolchain then decodes its labels but finds none by name: label_exists("562") is false on the published mini RefSeq index although decode() returns "562". Read the stored record whole and rebuild the buckets from the labels. From d5d26baa metagraph/src/annotation/representation/base/annotation.cpp. Co-authored-by: Gunnar Rätsch Co-Authored-By: Claude Opus 5.5 --- metagraph/src/annotation/representation/base/annotation.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/metagraph/src/annotation/representation/base/annotation.cpp b/metagraph/src/annotation/representation/base/annotation.cpp index 4cae5c7108..93083e5b05 100644 --- a/metagraph/src/annotation/representation/base/annotation.cpp +++ b/metagraph/src/annotation/representation/base/annotation.cpp @@ -54,7 +54,11 @@ bool LabelEncoder::load(std::istream &instream) { std::string version(7, '\0'); if (instream.read(version.data(), 7) && version == "LE-v2.0") { Deserializer deserializer(instream); - encode_label_ = VectorSet::deserialize(deserializer, true); + // Read the saved buckets to consume the entire record, then rebuild them: + // std::hash values can differ between the writer's and reader's toolchains. + auto labels = to_vector(VectorSet::deserialize(deserializer, true)); + encode_label_ = VectorSet(std::make_move_iterator(labels.begin()), + std::make_move_iterator(labels.end())); return instream.good(); } // backward compatibility From 0e7fc6e51514c08e0e46ee421e692c43d809778a Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sun, 1 Nov 2026 23:28:45 +0100 Subject: [PATCH 04/52] Follow every continuation of a label with an unlimited branch allowance MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit max_label_branches "unlimited" excluded every label on two or more successors, the opposite of what it means: an unlimited allowance never excludes a label. Found by the parity check against the prototype (ndm1_unlimited: 10 left paths instead of 29). Co-authored-by: Gunnar Rätsch Co-Authored-By: Claude Opus 5.5 --- metagraph/src/graph/traversal/walker.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/metagraph/src/graph/traversal/walker.cpp b/metagraph/src/graph/traversal/walker.cpp index e9afb7eda4..138a8997f5 100644 --- a/metagraph/src/graph/traversal/walker.cpp +++ b/metagraph/src/graph/traversal/walker.cpp @@ -815,8 +815,8 @@ std::optional Walker::process_item(ArmState &arm, Item &item, for (const Entry &src : item.state) { auto it = continuations.find(src.label); if (it != continuations.end() && it->second >= 2 - && (strategy_.max_label_branches == Strategy::kUnlimited - || src.branches + 1 > strategy_.max_label_branches)) { + && strategy_.max_label_branches != Strategy::kUnlimited + && src.branches + 1 > strategy_.max_label_branches) { excluded.insert(src.label); } } From 1483d511e422ca95f7304d93dce3b4daadd00920 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Mon, 2 Nov 2026 02:23:42 +0100 Subject: [PATCH 05/52] Keep g++-12 from a false array-bounds error in the walker tests GCC 12 (-Warray-bounds, -Werror) reports assigning a one-element initializer list to a std::vector as reading past the list; assign(1, x) says the same without the temporary array. Co-Authored-By: Claude Opus 5.5 --- metagraph/tests/graph/traversal/test_walker.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/metagraph/tests/graph/traversal/test_walker.cpp b/metagraph/tests/graph/traversal/test_walker.cpp index 9bfd59fab7..69a440e91f 100644 --- a/metagraph/tests/graph/traversal/test_walker.cpp +++ b/metagraph/tests/graph/traversal/test_walker.cpp @@ -952,7 +952,7 @@ TEST(WalkerDerive, IntersectionStartsFromTheCheapestRowOfTheFirstWindow) { for (size_t i = 0; i < labels.size(); ++i) { Column c = static_cast(encoder.encode(labels[i])); headers[c] = { "h_" + labels[i] }; - num_kmers[c] = { seqs[i].size() - kK + 1 }; + num_kmers[c].assign(1, seqs[i].size() - kK + 1); } annot::CoordToHeader cth(std::move(headers), std::move(num_kmers)); @@ -982,9 +982,9 @@ TEST(WalkerDerive, AmbiguousDerivedHeaderNamesAreRefused) { std::vector> headers(2); std::vector> num_kmers(2); headers[static_cast(encoder.encode("colA"))] = { "ACC1" }; - num_kmers[static_cast(encoder.encode("colA"))] = { (S + T).size() - kK + 1 }; + num_kmers[static_cast(encoder.encode("colA"))].assign(1, (S + T).size() - kK + 1); headers[static_cast(encoder.encode("colB"))] = { "ACC1" }; - num_kmers[static_cast(encoder.encode("colB"))] = { S.size() - kK + 1 }; + num_kmers[static_cast(encoder.encode("colB"))].assign(1, S.size() - kK + 1); annot::CoordToHeader cth(std::move(headers), std::move(num_kmers)); LabelOracle oracle(*anno, &cth); @@ -1021,9 +1021,9 @@ TEST(WalkerDerive, DerivedHeaderMustNotBeAColumnName) { std::vector> headers(2); std::vector> num_kmers(2); headers[enc.encode("ACC")] = { "decoy" }; - num_kmers[enc.encode("ACC")] = { in_acc.size() - kK + 1 }; + num_kmers[enc.encode("ACC")].assign(1, in_acc.size() - kK + 1); headers[enc.encode("B")] = { "ACC" }; - num_kmers[enc.encode("B")] = { in_b.size() - kK + 1 }; + num_kmers[enc.encode("B")].assign(1, in_b.size() - kK + 1); annot::CoordToHeader cth(std::move(headers), std::move(num_kmers)); LabelOracle oracle(*anno, &cth); From 6961a5953ced7190c1258fca5c93e64b9a5e38a4 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sun, 8 Nov 2026 20:30:22 +0100 Subject: [PATCH 06/52] Keep seed rank order and batch explicit label reads in /resolve Merging overlapping picks sorted the seeds by position. MAX_SUPPORT picks never overlap, so under the default policy merging only lost the best-first order. Merge LONGEST_FIRST picks only, give a merged pick the rank of its best member, and keep the seeds in rank order. Explicit labels decoded every row of the query in one call. Read them kResolveBatchRows at a time, as discovery does. Reject discover_max_labels = 0, and do not count cache hits for rows that are read again after a cache clear. Co-Authored-By: Claude Opus 5.5 --- .../src/graph/traversal/label_oracle.cpp | 5 +++- metagraph/src/graph/traversal/resolve.cpp | 29 ++++++++++++++----- .../tests/graph/traversal/test_resolve.cpp | 25 ++++++++++++++++ 3 files changed, 50 insertions(+), 9 deletions(-) diff --git a/metagraph/src/graph/traversal/label_oracle.cpp b/metagraph/src/graph/traversal/label_oracle.cpp index 60f1c89236..80de891ed3 100644 --- a/metagraph/src/graph/traversal/label_oracle.cpp +++ b/metagraph/src/graph/traversal/label_oracle.cpp @@ -246,18 +246,20 @@ std::vector LabelQuery::fetch(const std::vector missing; tsl::hopscotch_set requested; + size_t cache_hits = 0; for (node_index key : keys) { if (key == npos) continue; counters.rows_requested++; if (cache_.count(key)) { - counters.cache_hits++; + cache_hits++; } else if (requested.insert(key).second) { missing.push_back(key); } } // a clear in fetch_uncached may drop keys requested here: read them all again then if (cache_.size() + missing.size() > max_cache_size_) { + cache_hits = 0; missing.clear(); requested.clear(); for (node_index key : keys) { @@ -265,6 +267,7 @@ std::vector LabelQuery::fetch(const std::vector result(keys.size()); diff --git a/metagraph/src/graph/traversal/resolve.cpp b/metagraph/src/graph/traversal/resolve.cpp index e1bd6bf25f..4e3311952c 100644 --- a/metagraph/src/graph/traversal/resolve.cpp +++ b/metagraph/src/graph/traversal/resolve.cpp @@ -162,6 +162,8 @@ SupportProfile resolve_support(const LabelOracle &oracle, throw std::invalid_argument("Specify either an explicit label list or discovery, not both"); if (query.size() < profile.k) throw std::invalid_argument("Query shorter than k"); + if (options.discover && !options.discover_max_labels) + throw std::invalid_argument("Discovery needs discover_max_labels of at least 1"); profile.num_kmers = query.size() - profile.k + 1; const std::vector keys = oracle.keys_of_sequence(query); @@ -184,11 +186,16 @@ SupportProfile resolve_support(const LabelOracle &oracle, refs.push_back(oracle.resolve_label(name)); } LabelQuery query_labels(oracle, refs); - const auto hits = query_labels.fetch(keys); support.resize(refs.size()); - for (uint64_t i = 0; i < hits.size(); ++i) { - for (LabelId l : hits[i]) { - support[l].add(i); + // decode the rows kResolveBatchRows at a time, as discovery does + for (size_t begin = 0; begin < keys.size(); begin += kResolveBatchRows) { + const size_t end = std::min(keys.size(), begin + kResolveBatchRows); + const auto hits = query_labels.fetch( + std::vector(keys.begin() + begin, keys.begin() + end)); + for (size_t i = 0; i < hits.size(); ++i) { + for (LabelId l : hits[i]) { + support[l].add(begin + i); + } } } } else { @@ -322,6 +329,7 @@ std::vector covering_labels(const SupportProfile &profile, const KmerIn struct Picked { KmerInterval kmers; std::vector labels; + size_t rank; // the order in which the policy picked it }; // MAX_SUPPORT: repeatedly pick the interval of at least |min_kmers| k-mers contained in the @@ -419,7 +427,7 @@ std::vector pick_max_support(const SupportProfile &profile, uint64_t min } if (!best) break; - picked.push_back({ *best, covering_labels(profile, *best) }); + picked.push_back({ *best, covering_labels(profile, *best), picked.size() }); assert(picked.back().labels.size() == best_count); taken.push_back(*best); std::sort(taken.begin(), taken.end()); @@ -448,15 +456,17 @@ SeedSelection select_seeds(const SupportProfile &profile, continue; selection.num_eligible++; if (picked.size() < policy.max_seeds) - picked.push_back({ cand.kmers, cand.labels }); + picked.push_back({ cand.kmers, cand.labels, picked.size() }); } } else { picked = pick_max_support(profile, min_kmers, policy.max_seeds); selection.num_eligible = picked.size(); } - // merge overlapping picks onto their common interval - if (policy.merge_overlapping) { + // Merge overlapping picks onto their common interval, carried by the labels of both. + // Only LONGEST_FIRST picks can overlap: MAX_SUPPORT picks outside the earlier ones. + // A merged pick keeps the rank of its best member, and the picks stay in rank order. + if (policy.merge_overlapping && policy.policy == SelectionPolicy::LONGEST_FIRST) { std::sort(picked.begin(), picked.end(), [](const Picked &a, const Picked &b) { return a.kmers < b.kmers; }); std::vector merged; @@ -468,11 +478,14 @@ SeedSelection select_seeds(const SupportProfile &profile, if (common.size() >= min_kmers) { m.kmers = common; m.labels = covering_labels(profile, common); + m.rank = std::min(m.rank, p.rank); continue; } } merged.push_back(std::move(p)); } + std::sort(merged.begin(), merged.end(), + [](const Picked &a, const Picked &b) { return a.rank < b.rank; }); picked.swap(merged); } diff --git a/metagraph/tests/graph/traversal/test_resolve.cpp b/metagraph/tests/graph/traversal/test_resolve.cpp index 82a235ec53..f52b8ef104 100644 --- a/metagraph/tests/graph/traversal/test_resolve.cpp +++ b/metagraph/tests/graph/traversal/test_resolve.cpp @@ -128,6 +128,31 @@ TYPED_TEST(ResolveTest, ThreeStatesAndRuns) { EXPECT_EQ((std::vector{ "A", "B" }), sel.seeds[0].labels); EXPECT_EQ(2u, sel.seeds[0].population.supporting_total); + // longest_first with merging: [95,110){A} and [100,190){B} merge onto + // [100,110){A,B}, which keeps the rank of [100,190) and so comes first + pol = SelectionPolicy(); + pol.policy = SelectionPolicy::LONGEST_FIRST; + pol.max_seeds = 3; + sel = select_seeds(profile, q, pol, canonical); + ASSERT_EQ(2u, sel.seeds.size()); + EXPECT_EQ((KmerInterval{ 100, 110 }), sel.seeds[0].kmers); + EXPECT_EQ((std::vector{ "A", "B" }), sel.seeds[0].labels); + EXPECT_EQ((KmerInterval{ 0, 70 }), sel.seeds[1].kmers); + EXPECT_EQ((std::vector{ "A" }), sel.seeds[1].labels); + + // max_support picks never overlap, so merging leaves them and their order as they are + pol = SelectionPolicy(); + pol.max_seeds = 3; + pol.merge_overlapping = false; + auto unmerged = select_seeds(profile, q, pol, canonical); + pol.merge_overlapping = true; + sel = select_seeds(profile, q, pol, canonical); + ASSERT_EQ(unmerged.seeds.size(), sel.seeds.size()); + for (size_t i = 0; i < sel.seeds.size(); ++i) { + EXPECT_EQ(unmerged.seeds[i].kmers, sel.seeds[i].kmers) << i; + EXPECT_EQ(unmerged.seeds[i].labels, sel.seeds[i].labels) << i; + } + // max_labels_per_seed drops deterministically and reports the drop pol.max_labels_per_seed = 1; sel = select_seeds(profile, q, pol, canonical); From 546ea5c644484c2408962b63de81f55353f92c35 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sun, 8 Nov 2026 21:11:36 +0100 Subject: [PATCH 07/52] Refuse ambiguous label names A name that is both an annotation column and a sequence header resolved to the column, and a header of several records to its first record, without saying so. Refuse both with std::invalid_argument (a 400 on the server). CoordToHeader marks the headers of several records, can build its reverse index ahead of time (index_headers), and answers lookups outside its lock once the index is built. Co-Authored-By: Claude Opus 5.5 --- metagraph/src/annotation/coord_to_header.cpp | 27 +++++++++++++++---- metagraph/src/annotation/coord_to_header.hpp | 12 +++++++-- .../src/graph/traversal/label_oracle.cpp | 14 ++++++++-- .../src/graph/traversal/label_oracle.hpp | 5 ++-- .../graph/traversal/test_label_oracle.cpp | 17 +++++++++++- 5 files changed, 63 insertions(+), 12 deletions(-) diff --git a/metagraph/src/annotation/coord_to_header.cpp b/metagraph/src/annotation/coord_to_header.cpp index e03ddec80f..5637c5e1bb 100644 --- a/metagraph/src/annotation/coord_to_header.cpp +++ b/metagraph/src/annotation/coord_to_header.cpp @@ -39,13 +39,13 @@ struct HeaderOfRecordEqual { // header -> (column, seq_id) of its first occurrence in column order struct CoordToHeader::HeaderIndex { tsl::hopscotch_map> map; + tsl::hopscotch_set repeated; // the headers of several records }; CoordToHeader::CoordToHeader() {} CoordToHeader::~CoordToHeader() {} -std::optional> -CoordToHeader::find_header(std::string_view header) const { +const CoordToHeader::HeaderIndex& CoordToHeader::header_index() const { std::lock_guard lock(header_index_mutex_); if (!header_index_) { auto index = std::make_unique(); @@ -56,17 +56,34 @@ CoordToHeader::find_header(std::string_view header) const { index->map.reserve(total); for (Column col = 0; col < headers_.size(); ++col) { for (size_t s = 0; s < headers_[col].size(); ++s) { - index->map.try_emplace(std::string_view(headers_[col][s]), col, s); + const std::string_view header(headers_[col][s]); + if (!index->map.try_emplace(header, col, s).second) + index->repeated.insert(header); } } header_index_ = std::move(index); } - auto it = header_index_->map.find(header); - if (it == header_index_->map.end()) + // the index is immutable once built, so lookups run outside the lock + return *header_index_; +} + +void CoordToHeader::index_headers() const { + header_index(); +} + +std::optional> +CoordToHeader::find_header(std::string_view header) const { + const HeaderIndex &index = header_index(); + auto it = index.map.find(header); + if (it == index.map.end()) return std::nullopt; return it->second; } +bool CoordToHeader::is_repeated_header(std::string_view header) const { + return header_index().repeated.count(header); +} + void CoordToHeader::find_duplicate_headers(const std::vector &headers, const std::function &callback) { tsl::hopscotch_set diff --git a/metagraph/src/annotation/coord_to_header.hpp b/metagraph/src/annotation/coord_to_header.hpp index 6fe0f72651..2b9a164a89 100644 --- a/metagraph/src/annotation/coord_to_header.hpp +++ b/metagraph/src/annotation/coord_to_header.hpp @@ -91,11 +91,18 @@ class CoordToHeader { /** * The (column, seq_id) of a FASTA header: its first occurrence in column order (a * header is unique within a column, not within the index), or nothing. The reverse - * index is built on the first call and owned by this object, so it lives and dies with - * the headers it indexes. load() drops it. Thread-safe. + * index is built by index_headers() or on the first lookup, and owned by this object, + * so it lives and dies with the headers it indexes. load() drops it. Thread-safe, + * except against a concurrent load(). */ std::optional> find_header(std::string_view header) const; + // Whether several records carry |header|, in one column or in several. + bool is_repeated_header(std::string_view header) const; + + // Builds the reverse header index now, so that no lookup pays for it. + void index_headers() const; + /** * The records of |headers| whose header repeats an earlier record's, in record order: * |callback|(first, repeat) with the 0-based indices of the first occurrence and of the @@ -110,6 +117,7 @@ class CoordToHeader { private: struct HeaderIndex; + const HeaderIndex& header_index() const; std::vector> headers_; std::vector coord_offsets_; diff --git a/metagraph/src/graph/traversal/label_oracle.cpp b/metagraph/src/graph/traversal/label_oracle.cpp index 80de891ed3..e300fe2481 100644 --- a/metagraph/src/graph/traversal/label_oracle.cpp +++ b/metagraph/src/graph/traversal/label_oracle.cpp @@ -125,14 +125,24 @@ std::pair LabelOracle::sequence_coords(Column column, uint64 } LabelRef LabelOracle::resolve_label(const std::string &name) const { - if (auto column = find_column(name)) { + auto column = find_column(name); + auto header = find_header(name); + if (column && header) { + throw std::invalid_argument("Label '" + name + "' is ambiguous: it names both an " + "annotation column and a sequence header"); + } + if (header && coord_to_header_->is_repeated_header(name)) { + throw std::invalid_argument("Label '" + name + "' is ambiguous: several sequence " + "records have this header"); + } + if (column) { LabelRef ref; ref.kind = LabelKind::COLUMN; ref.column = *column; ref.name = name; return ref; } - if (auto header = find_header(name)) + if (header) return *header; throw std::invalid_argument("Label not found in the annotation: '" + name + "'" + (coord_to_header_ ? "" : " (no sequence header index loaded)")); diff --git a/metagraph/src/graph/traversal/label_oracle.hpp b/metagraph/src/graph/traversal/label_oracle.hpp index a92cf5e224..1ebea45d91 100644 --- a/metagraph/src/graph/traversal/label_oracle.hpp +++ b/metagraph/src/graph/traversal/label_oracle.hpp @@ -81,8 +81,9 @@ class LabelOracle { std::pair map_coord(Column column, uint64_t coord) const; // The first and last column coordinates of sequence |seq_id|. Requires a CoordToHeader. std::pair sequence_coords(Column column, uint64_t seq_id) const; - // A column name first, then a header. Throws std::invalid_argument naming the label if - // neither exists. + // A column name or a header. Throws std::invalid_argument naming the label if neither + // exists, or if it is ambiguous: both a column and a header, or the header of several + // records. LabelRef resolve_label(const std::string &name) const; // Batched row access (rows = annotation key - 1). Updates the counters. diff --git a/metagraph/tests/graph/traversal/test_label_oracle.cpp b/metagraph/tests/graph/traversal/test_label_oracle.cpp index a3c75d8054..123d4551af 100644 --- a/metagraph/tests/graph/traversal/test_label_oracle.cpp +++ b/metagraph/tests/graph/traversal/test_label_oracle.cpp @@ -334,6 +334,17 @@ TYPED_TEST(LabelOracleCoordTest, HeaderLabels) { } EXPECT_GT(oracle.counters().tuple_rows_fetched, 0u); + // a name of both a column and a header, or the header of several records, is ambiguous + std::vector> clash(2); + std::vector> clash_kmers(2); + clash[col_f] = { "G", "acc2" }; clash_kmers[col_f] = { n1, n2 }; + clash[col_g] = { "acc2" }; clash_kmers[col_g] = { n3 }; + annot::CoordToHeader clash_cth(std::move(clash), std::move(clash_kmers)); + LabelOracle clash_oracle(*anno, &clash_cth); + EXPECT_THROW(clash_oracle.resolve_label("G"), std::invalid_argument); + EXPECT_THROW(clash_oracle.resolve_label("acc2"), std::invalid_argument); + EXPECT_EQ(LabelKind::COLUMN, clash_oracle.resolve_label("F").kind); + // header labels need a CoordToHeader LabelOracle no_headers(*anno); EXPECT_THROW(no_headers.resolve_label("acc1"), std::invalid_argument); @@ -359,10 +370,14 @@ TEST(CoordToHeaderIndex, FindHeader) { EXPECT_FALSE(oracle.find_header("W")); } - // a header repeated in a second column resolves to its first column + // a header repeated in a second column resolves to its first column, and is marked annot::CoordToHeader two({ { "A", "B" }, { "B", "C" } }, { { 1, 1 }, { 1, 1 } }); + two.index_headers(); EXPECT_EQ(std::make_pair(Column(0), size_t(1)), *two.find_header("B")); EXPECT_EQ(std::make_pair(Column(1), size_t(1)), *two.find_header("C")); + EXPECT_TRUE(two.is_repeated_header("B")); + EXPECT_FALSE(two.is_repeated_header("C")); + EXPECT_FALSE(two.is_repeated_header("W")); // load() replaces the headers and the index with them const std::string base = test_dump_dir() + "/coord_to_header_index"; From efa09d8b605cacfac9b80e200f4d9958322ac2ff Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sun, 8 Nov 2026 21:11:53 +0100 Subject: [PATCH 08/52] Check derived header labels without resolving them resolve_label now refuses an ambiguous name. Derivation checks the same two conditions directly, so an ambiguous derived header still ends the seed with ambiguous_header instead of failing the request. Co-Authored-By: Claude Opus 5.5 --- metagraph/src/graph/traversal/walker.cpp | 14 ++++++-------- 1 file changed, 6 insertions(+), 8 deletions(-) diff --git a/metagraph/src/graph/traversal/walker.cpp b/metagraph/src/graph/traversal/walker.cpp index 138a8997f5..07e9cb4c12 100644 --- a/metagraph/src/graph/traversal/walker.cpp +++ b/metagraph/src/graph/traversal/walker.cpp @@ -476,19 +476,17 @@ void Walker::derive_seed_labels(const std::vector &keys, : oracle_.header_name(key.first, key.second); }; if (kind == LabelKind::HEADER) { - // A derived list must be resubmittable verbatim as an explicit one, which resolves a - // name to a column first, then to the first column holding the header + // A derived list must be resubmittable verbatim as an explicit one, and + // resolve_label refuses a header that also names a column or another record for (const Key &key : live) { const std::string &name = name_of(key); - const LabelRef back = oracle_.resolve_label(name); // the header exists: no throw - if (back.kind != LabelKind::HEADER || back.column != key.first - || back.seq_id != key.second) { + const bool column_too = oracle_.find_column(name).has_value(); + if (column_too || oracle_.coord_to_header()->is_repeated_header(name)) { throw SeedDerivationError(SeedDerivationError::AMBIGUOUS_HEADER, "The labels derived from the seed are not resubmittable: the sequence " "header '" + name + "' " - + (back.kind == LabelKind::COLUMN - ? "is also the name of an annotation column" - : "also occurs in another annotation column") + + (column_too ? "is also the name of an annotation column" + : "is also the header of another sequence record") + ". Set seed_label_kind to \"column\" or name the labels explicitly.", 0, 0, name); } From 7547cf5daf541a4959e16a9888fbb63baa8787c0 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sun, 8 Nov 2026 21:19:53 +0100 Subject: [PATCH 09/52] Bound discovery and skip columns without sequence records Discovery kept every label it saw before keeping the best max_labels, so its memory grew with labels and runs. Refuse a discovery that sees more than 100,000 distinct labels (kMaxDiscoveredLabels). Header discovery mapped the coordinates of every column, and a column without sequence records threw out_of_range, failing the request. Skip such columns: they have no header labels (LabelOracle::has_header_records). Co-Authored-By: Claude Opus 5.5 --- metagraph/src/graph/traversal/label_oracle.cpp | 5 +++++ metagraph/src/graph/traversal/label_oracle.hpp | 3 +++ metagraph/src/graph/traversal/resolve.cpp | 7 +++++++ metagraph/src/graph/traversal/resolve.hpp | 2 ++ .../graph/traversal/test_label_oracle.cpp | 18 ++++++++++++++++++ 5 files changed, 35 insertions(+) diff --git a/metagraph/src/graph/traversal/label_oracle.cpp b/metagraph/src/graph/traversal/label_oracle.cpp index e300fe2481..98c525355f 100644 --- a/metagraph/src/graph/traversal/label_oracle.cpp +++ b/metagraph/src/graph/traversal/label_oracle.cpp @@ -111,6 +111,11 @@ const std::string& LabelOracle::header_name(Column column, uint64_t seq_id) cons return coord_to_header_->get_headers(column).at(seq_id); } +bool LabelOracle::has_header_records(Column column) const { + return coord_to_header_ && column < coord_to_header_->num_columns() + && coord_to_header_->num_sequences(column); +} + std::pair LabelOracle::map_coord(Column column, uint64_t coord) const { assert(coord_to_header_); counters_.coords_mapped++; diff --git a/metagraph/src/graph/traversal/label_oracle.hpp b/metagraph/src/graph/traversal/label_oracle.hpp index 1ebea45d91..29c90bb8c9 100644 --- a/metagraph/src/graph/traversal/label_oracle.hpp +++ b/metagraph/src/graph/traversal/label_oracle.hpp @@ -77,6 +77,9 @@ class LabelOracle { // The first (column, sequence) holding a FASTA header (CoordToHeader::find_header) std::optional find_header(const std::string &name) const; const std::string& header_name(Column column, uint64_t seq_id) const; + // Whether |column| has sequence records, and so header labels (false without a + // CoordToHeader). + bool has_header_records(Column column) const; // (seq_id, local coordinate) of a column coordinate. Requires a CoordToHeader. std::pair map_coord(Column column, uint64_t coord) const; // The first and last column coordinates of sequence |seq_id|. Requires a CoordToHeader. diff --git a/metagraph/src/graph/traversal/resolve.cpp b/metagraph/src/graph/traversal/resolve.cpp index 4e3311952c..ed7d62c1a8 100644 --- a/metagraph/src/graph/traversal/resolve.cpp +++ b/metagraph/src/graph/traversal/resolve.cpp @@ -121,6 +121,11 @@ void discover_labels(const LabelOracle &oracle, auto of = [&](Column c, uint64_t seq_id) -> LabelSupport& { auto [it, inserted] = index.try_emplace({ c, seq_id }, found->size()); if (inserted) { + if (found->size() == kMaxDiscoveredLabels) { + throw std::invalid_argument("Discovery found more than " + + std::to_string(kMaxDiscoveredLabels) + " distinct labels. " + "Use a shorter query, or name the labels explicitly."); + } found->emplace_back(); ids->emplace_back(c, seq_id); } @@ -140,6 +145,8 @@ void discover_labels(const LabelOracle &oracle, const auto rows = oracle.get_row_tuples(anno_rows(batch)); for (size_t j = 0; j < rows.size(); ++j) { for (const auto &[c, coords] : rows[j]) { + if (!oracle.has_header_records(c)) + continue; // no records, so no header labels for (uint64_t coord : coords) { of(c, oracle.map_coord(c, coord).first).add(pos[begin + j]); } diff --git a/metagraph/src/graph/traversal/resolve.hpp b/metagraph/src/graph/traversal/resolve.hpp index 9188997dc6..e8e2863645 100644 --- a/metagraph/src/graph/traversal/resolve.hpp +++ b/metagraph/src/graph/traversal/resolve.hpp @@ -29,6 +29,8 @@ struct KmerInterval { // The rows a discovery decodes at once constexpr size_t kResolveBatchRows = 4096; +// The most distinct labels a discovery may see before it keeps the best |discover_max_labels| +constexpr size_t kMaxDiscoveredLabels = 100'000; struct ResolveOptions { // exactly one of |labels| (explicit) or |discover| must be set diff --git a/metagraph/tests/graph/traversal/test_label_oracle.cpp b/metagraph/tests/graph/traversal/test_label_oracle.cpp index 123d4551af..bab843a1b1 100644 --- a/metagraph/tests/graph/traversal/test_label_oracle.cpp +++ b/metagraph/tests/graph/traversal/test_label_oracle.cpp @@ -8,6 +8,7 @@ #include "tests/annotation/test_annotated_dbg_helpers.hpp" #include "graph/traversal/label_oracle.hpp" +#include "graph/traversal/resolve.hpp" #include "graph/annotated_dbg.hpp" #include "graph/representation/canonical_dbg.hpp" #include "graph/representation/hash/dbg_sshash.hpp" @@ -345,6 +346,23 @@ TYPED_TEST(LabelOracleCoordTest, HeaderLabels) { EXPECT_THROW(clash_oracle.resolve_label("acc2"), std::invalid_argument); EXPECT_EQ(LabelKind::COLUMN, clash_oracle.resolve_label("F").kind); + // a column without sequence records has no header labels: header discovery skips it + std::vector> partial(2); + std::vector> partial_kmers(2); + partial[col_f] = { "acc1", "acc2" }; partial_kmers[col_f] = { n1, n2 }; + annot::CoordToHeader partial_cth(std::move(partial), std::move(partial_kmers)); + LabelOracle partial_oracle(*anno, &partial_cth); + EXPECT_TRUE(partial_oracle.has_header_records(col_f)); + EXPECT_FALSE(partial_oracle.has_header_records(col_g)); + ResolveOptions discover; + discover.discover = true; + discover.discover_kind = LabelKind::HEADER; + auto profile = resolve_support(partial_oracle, s3, discover); + ASSERT_FALSE(profile.labels.empty()); + for (const auto &support : profile.labels) { + EXPECT_EQ(col_f, support.label.column); + } + // header labels need a CoordToHeader LabelOracle no_headers(*anno); EXPECT_THROW(no_headers.resolve_label("acc1"), std::invalid_argument); From b02de7e81384e479aa36f18e33712ce530679661 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sun, 8 Nov 2026 21:20:11 +0100 Subject: [PATCH 10/52] Skip columns without sequence records in header derivation A column without sequence records threw out_of_range when the walker derived header labels from a seed. Skip it, as discovery does. The derivation test's explicit request for a header of two records now expects the refusal of an ambiguous name. Co-Authored-By: Claude Opus 5.5 --- metagraph/src/graph/traversal/walker.cpp | 2 ++ metagraph/tests/graph/traversal/test_walker.cpp | 4 ++-- 2 files changed, 4 insertions(+), 2 deletions(-) diff --git a/metagraph/src/graph/traversal/walker.cpp b/metagraph/src/graph/traversal/walker.cpp index 07e9cb4c12..b08163c024 100644 --- a/metagraph/src/graph/traversal/walker.cpp +++ b/metagraph/src/graph/traversal/walker.cpp @@ -432,6 +432,8 @@ void Walker::derive_seed_labels(const std::vector &keys, })) { continue; } + if (!oracle_.has_header_records(c)) + continue; // no records, so no header labels for (uint64_t coord : coords) { cur.emplace_back(c, oracle_.map_coord(c, coord).first); } diff --git a/metagraph/tests/graph/traversal/test_walker.cpp b/metagraph/tests/graph/traversal/test_walker.cpp index 69a440e91f..f9e36a1bb0 100644 --- a/metagraph/tests/graph/traversal/test_walker.cpp +++ b/metagraph/tests/graph/traversal/test_walker.cpp @@ -1005,9 +1005,9 @@ TEST(WalkerDerive, AmbiguousDerivedHeaderNamesAreRefused) { ASSERT_EQ(2u, by_column.label_dict.size()); EXPECT_EQ("colA", by_column.label_dict[0].name); EXPECT_EQ("colB", by_column.label_dict[1].name); + // named explicitly, the header of two records is refused as ambiguous seed.labels = { "ACC1" }; - auto one = traverse_seed(oracle, seed, Strategy()); - EXPECT_EQ(1u, one.label_dict.size()); + EXPECT_THROW(traverse_seed(oracle, seed, Strategy()), std::invalid_argument); } // a derived header that an explicit list would resolve to a column is refused From 74b4c4644d4de1ac8c8e144ecb2700f069a13ee1 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sun, 8 Nov 2026 21:27:10 +0100 Subject: [PATCH 11/52] Keep the label cache bounded and make CoordToHeader::load all-or-nothing LabelQuery cleared an over-full cache and then stored the whole batch in it, so a batch with more keys than max_cache_size grew the cache past its bound. Answer such a batch from a temporary map. CoordToHeader::load replaced its vectors in place, so a failed load left a half-loaded object. Read into new vectors and swap them in, dropping the header index in the same step. Co-Authored-By: Claude Opus 5.5 --- metagraph/src/annotation/coord_to_header.cpp | 20 +++++++-------- .../src/graph/traversal/label_oracle.cpp | 25 +++++++++++-------- .../src/graph/traversal/label_oracle.hpp | 6 +++-- 3 files changed, 29 insertions(+), 22 deletions(-) diff --git a/metagraph/src/annotation/coord_to_header.cpp b/metagraph/src/annotation/coord_to_header.cpp index 5637c5e1bb..63a37fddfa 100644 --- a/metagraph/src/annotation/coord_to_header.cpp +++ b/metagraph/src/annotation/coord_to_header.cpp @@ -191,21 +191,21 @@ bool CoordToHeader::load(const std::string &filename_base) { utils::file_read_failure_detail(path)); return false; } - { - // the index views the headers about to be replaced - std::lock_guard lock(header_index_mutex_); - header_index_.reset(); - } - try { + // read into new vectors, so that a failed load leaves this object as it was uint64_t num_columns = load_number(*in); - headers_.resize(num_columns); - coord_offsets_.resize(num_columns); + std::vector> headers(num_columns); + std::vector coord_offsets(num_columns); for (uint64_t i = 0; i < num_columns; ++i) { - load_string_vector(*in, &headers_[i]); - coord_offsets_[i].load(*in); + load_string_vector(*in, &headers[i]); + coord_offsets[i].load(*in); } + // the index views the headers being replaced: drop it with them + std::lock_guard lock(header_index_mutex_); + header_index_.reset(); + headers_.swap(headers); + coord_offsets_.swap(coord_offsets); return true; } catch (const std::exception &e) { diff --git a/metagraph/src/graph/traversal/label_oracle.cpp b/metagraph/src/graph/traversal/label_oracle.cpp index 98c525355f..d8796c3d9b 100644 --- a/metagraph/src/graph/traversal/label_oracle.cpp +++ b/metagraph/src/graph/traversal/label_oracle.cpp @@ -70,7 +70,9 @@ std::vector LabelOracle::keys_of_sequence(std::string_view sequence) return keys; keys.reserve(sequence.size() - get_k() + 1); // map_to_nodes returns annotation keys in every regime: the node itself (BASIC), - // the primary node (CanonicalDBG), the canonical node (native CANONICAL graphs) + // the primary node (CanonicalDBG), the canonical node (native CANONICAL graphs). + // Every graph calls them in query order, although the base class promises that only + // for map_to_nodes_sequentially; RegimeAndKeyMapping checks it per k-mer. graph_->map_to_nodes(sequence, [&](node_index key) { keys.push_back(key); }); counters_.keys_mapped += keys.size(); return keys; @@ -233,12 +235,10 @@ void LabelQuery::hits_from_tuples(const MultiIntMatrix::RowTuples &row, NodeHits hits->erase(std::unique(hits->begin(), hits->end()), hits->end()); } -void LabelQuery::fetch_uncached(const std::vector &keys) { +void LabelQuery::read_hits(const std::vector &keys, + tsl::hopscotch_map *out) const { if (keys.empty()) return; - if (cache_.size() + keys.size() > max_cache_size_) - cache_.clear(); - std::vector rows; rows.reserve(keys.size()); for (node_index key : keys) { @@ -247,12 +247,12 @@ void LabelQuery::fetch_uncached(const std::vector &keys) { if (tuples_) { auto tuple_rows = oracle_.get_row_tuples(rows); for (size_t i = 0; i < keys.size(); ++i) { - hits_from_tuples(tuple_rows[i], &cache_[keys[i]]); + hits_from_tuples(tuple_rows[i], &(*out)[keys[i]]); } } else { auto set_rows = oracle_.get_rows(rows); for (size_t i = 0; i < keys.size(); ++i) { - hits_from_row(set_rows[i], &cache_[keys[i]]); + hits_from_row(set_rows[i], &(*out)[keys[i]]); } } } @@ -272,9 +272,10 @@ std::vector LabelQuery::fetch(const std::vector max_cache_size_) { cache_hits = 0; + cache_.clear(); missing.clear(); requested.clear(); for (node_index key : keys) { @@ -283,12 +284,16 @@ std::vector LabelQuery::fetch(const std::vector uncached; + auto *found = missing.size() > max_cache_size_ ? &uncached : &cache_; + read_hits(missing, found); std::vector result(keys.size()); for (size_t i = 0; i < keys.size(); ++i) { if (keys[i] != npos) - result[i] = cache_.at(keys[i]); + result[i] = found->at(keys[i]); } return result; } diff --git a/metagraph/src/graph/traversal/label_oracle.hpp b/metagraph/src/graph/traversal/label_oracle.hpp index 29c90bb8c9..d57bb1bd17 100644 --- a/metagraph/src/graph/traversal/label_oracle.hpp +++ b/metagraph/src/graph/traversal/label_oracle.hpp @@ -116,7 +116,8 @@ class LabelOracle { * Batched membership of a fixed set of permitted labels at annotation keys. Column * labels are read from whole rows; header labels (one sequence of a column) from tuple * rows, a coordinate counting for a label when it falls into that label's sequence. - * Results are cached per key, up to |max_cache_size| keys, then the cache is cleared. + * Results are cached per key, up to |max_cache_size| keys, then the cache is cleared. A + * batch with more keys than that is answered without caching it. */ class LabelQuery { public: @@ -143,7 +144,8 @@ class LabelQuery { LabelId label; }; - void fetch_uncached(const std::vector &keys); + void read_hits(const std::vector &keys, + tsl::hopscotch_map *out) const; void hits_from_row(const annot::matrix::BinaryMatrix::SetBitPositions &row, NodeHits *hits) const; void hits_from_tuples(const annot::matrix::MultiIntMatrix::RowTuples &row, From 32426c31346bcb9e2f40d72afd431f657719f186 Mon Sep 17 00:00:00 2001 From: Oleksandr Kulkov Date: Sun, 8 Nov 2026 22:25:25 +0100 Subject: [PATCH 12/52] Note that the label encoder's saved buckets can be dropped load() rebuilds the hash buckets from the labels (#725), so the buckets saved in each LE-v2.0 record are unused. Dropping them needs a format version that earlier builds cannot read, so it waits for the next change to this format. Co-Authored-By: Claude Opus 5.5 --- metagraph/src/annotation/representation/base/annotation.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/metagraph/src/annotation/representation/base/annotation.cpp b/metagraph/src/annotation/representation/base/annotation.cpp index 93083e5b05..8a5cb93a68 100644 --- a/metagraph/src/annotation/representation/base/annotation.cpp +++ b/metagraph/src/annotation/representation/base/annotation.cpp @@ -30,6 +30,8 @@ size_t LabelEncoder