support for multiple successors at forks

metagraph/src/annotation/annotation_converters.cpp

@@ -1157,7 +1157,8 @@ void convert_to_row_diff(const std::vector<std::string> &files,
 
     if (construction_stage == RowDiffStage::CONVERT) {
         assign_anchors(graph_fname, graph_fname, out_dir, max_path_length,
-                       ".row_reduction", get_num_threads());
+                       ".row_reduction", get_num_threads(),
+                       with_values || with_coordinates);
 
         const std::string anchors_fname = graph_fname + kRowDiffAnchorExt;
         if (!fs::exists(anchors_fname)) {
@@ -1172,6 +1173,20 @@ void convert_to_row_diff(const std::vector<std::string> &files,
             return;
         }
         mem_bytes -= anchor_size;
+
+        const std::string rd_succ_fname = graph_fname + kRowDiffForkSuccExt;
+        if (!fs::exists(rd_succ_fname)) {
+            logger->error("Can't find row-diff successor bitmap at {}", rd_succ_fname);
+            exit(1);
+        }
+        uint64_t rd_succ_size = fs::file_size(rd_succ_fname);
+        if (rd_succ_size > mem_bytes) {
+            logger->warn("row-diff successor bitmap ({} MiB) is larger than"
+                         " the memory allocated ({} MiB). Reserve more RAM.",
+                         rd_succ_size >> 20, mem_bytes >> 20);
+            return;
+        }
+        mem_bytes -= rd_succ_size;
     }
 
     if (!files.size())

metagraph/src/annotation/binary_matrix/row_diff/row_diff.cpp

@@ -37,13 +37,15 @@ void IRowDiff::load_fork_succ(const std::string &filename) {
     fork_succ_.load(f);
 }
 
-std::pair<std::vector<BinaryMatrix::Row>, std::vector<std::vector<size_t>>>
+std::pair<std::vector<BinaryMatrix::Row>, std::vector<std::vector<std::pair<size_t, size_t>>>>
 IRowDiff::get_rd_ids(const std::vector<BinaryMatrix::Row> &row_ids) const {
     assert(graph_ && "graph must be loaded");
-    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     using Row = BinaryMatrix::Row;
 
+    const size_t RD_PATH_RESERVE_SIZE = 2;
+
     // diff rows annotating nodes along the row-diff paths
     std::vector<Row> rd_ids;
     rd_ids.reserve(row_ids.size() * RD_PATH_RESERVE_SIZE);
@@ -55,38 +57,50 @@ IRowDiff::get_rd_ids(const std::vector<BinaryMatrix::Row> &row_ids) const {
     // Truncated row-diff paths, indexes to |rd_rows|.
     // The last index in each path points to an anchor or to a row which had
     // been reached before, and thus, will be reconstructed before this one.
-    std::vector<std::vector<size_t>> rd_paths_trunc(row_ids.size());
-
-    const graph::boss::BOSS &boss = graph_->get_boss();
-    const bit_vector &rd_succ = fork_succ_.size() ? fork_succ_ : boss.get_last();
+    std::vector<std::vector<std::pair<size_t, size_t>>> rd_paths_trunc(row_ids.size());
 
     for (size_t i = 0; i < row_ids.size(); ++i) {
-        Row row = row_ids[i];
-
-        graph::boss::BOSS::edge_index boss_edge = graph_->kmer_to_boss_index(
-                graph::AnnotatedSequenceGraph::anno_to_graph_index(row));
-
-        while (true) {
-            row = graph::AnnotatedSequenceGraph::graph_to_anno_index(
-                    graph_->boss_to_kmer_index(boss_edge));
-
+        std::vector<std::pair<size_t, size_t>> &rd_path = rd_paths_trunc[i];
+
+        std::vector<size_t> path;
+        Vector<std::pair<size_t, Row>> queue;
+        queue.emplace_back(0, row_ids[i]);
+
+        while (queue.size()) {
+            size_t depth = queue.back().first;
+            Row row = queue.back().second;
+            queue.pop_back();
+            while (depth < path.size()) {
+                assert(path.size() > 1);
+                rd_path.emplace_back(*(path.rbegin() + 1), *path.rbegin());
+                path.pop_back();
+            }
             auto [it, is_new] = node_to_rd.try_emplace(row, rd_ids.size());
-            rd_paths_trunc[i].push_back(it.value());
-
+            path.push_back(it.value());
             // If a node had been reached before, we interrupt the diff path.
             // The annotation for that node will have been reconstructed earlier
             // than for other nodes in this path as well. Thus, we will start
             // reconstruction from that node and don't need its successors.
             if (!is_new)
-                break;
+                continue;
 
             rd_ids.push_back(row);
 
             if (anchor_[row])
-                break;
+                continue;
+
+            auto node = graph::AnnotatedSequenceGraph::anno_to_graph_index(row);
+            graph_->call_row_diff_successors(node, fork_succ_, [&](auto succ) {
+                queue.emplace_back(depth + 1, graph::AnnotatedSequenceGraph::graph_to_anno_index(succ));
+            });
+        }
 
-            boss_edge = boss.row_diff_successor(boss_edge, rd_succ);
+        while (path.size() > 1) {
+            rd_path.emplace_back(*(path.rbegin() + 1), *path.rbegin());
+            path.pop_back();
         }
+        assert(path.size());
+        rd_path.emplace_back(-1, path[0]);
     }
 
     return std::make_pair(std::move(rd_ids), std::move(rd_paths_trunc));

metagraph/src/annotation/binary_matrix/row_diff/row_diff.hpp

@@ -13,7 +13,6 @@
 #include "common/logger.hpp"
 #include "common/utils/template_utils.hpp"
 #include "graph/annotated_dbg.hpp"
-#include "graph/representation/succinct/boss.hpp"
 #include "graph/representation/succinct/dbg_succinct.hpp"
 
 
@@ -44,7 +43,7 @@ class IRowDiff {
 
   protected:
     // get row-diff paths starting at |row_ids|
-    std::pair<std::vector<BinaryMatrix::Row>, std::vector<std::vector<size_t>>>
+    std::pair<std::vector<BinaryMatrix::Row>, std::vector<std::vector<std::pair<size_t, size_t>>>>
     get_rd_ids(const std::vector<BinaryMatrix::Row> &row_ids) const;
 
     const graph::DBGSuccinct *graph_ = nullptr;
@@ -115,7 +114,7 @@ template <class BaseMatrix>
 bool RowDiff<BaseMatrix>::get(Row row, Column column) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
-    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     SetBitPositions set_bits = get_row(row);
     SetBitPositions::iterator v = std::lower_bound(set_bits.begin(), set_bits.end(), column);
@@ -130,9 +129,9 @@ template <class BaseMatrix>
 std::vector<BinaryMatrix::Row> RowDiff<BaseMatrix>::get_column(Column column) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     const graph::boss::BOSS &boss = graph_->get_boss();
-    assert(!fork_succ_.size() || fork_succ_.size() == boss.get_last().size());
 
     std::vector<Row> result;
     // TODO: implement a more efficient algorithm
@@ -151,21 +150,16 @@ template <class BaseMatrix>
 BinaryMatrix::SetBitPositions RowDiff<BaseMatrix>::get_row(Row row) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
-    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     Vector<uint64_t> result = diffs_.get_row(row);
     std::sort(result.begin(), result.end());
 
-    uint64_t boss_edge = graph_->kmer_to_boss_index(
-            graph::AnnotatedSequenceGraph::anno_to_graph_index(row));
-    const graph::boss::BOSS &boss = graph_->get_boss();
-    const bit_vector &rd_succ = fork_succ_.size() ? fork_succ_ : boss.get_last();
+    auto node = graph::AnnotatedSequenceGraph::anno_to_graph_index(row);
 
     while (!anchor_[row]) {
-        boss_edge = boss.row_diff_successor(boss_edge, rd_succ);
-
-        row = graph::AnnotatedSequenceGraph::graph_to_anno_index(
-                graph_->boss_to_kmer_index(boss_edge));
+        node = graph_->row_diff_successor(node, fork_succ_);
+        row = graph::AnnotatedSequenceGraph::graph_to_anno_index(node);
 
         auto diff_row = diffs_.get_row(row);
         std::sort(diff_row.begin(), diff_row.end());
@@ -180,7 +174,7 @@ std::vector<BinaryMatrix::SetBitPositions>
 RowDiff<BaseMatrix>::get_rows(const std::vector<Row> &row_ids) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
-    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     // diff rows annotating nodes along the row-diff paths
     std::vector<Row> rd_ids;
@@ -199,19 +193,10 @@ RowDiff<BaseMatrix>::get_rows(const std::vector<Row> &row_ids) const {
     // been reached before, and thus, will be reconstructed before this one.
     std::vector<std::vector<size_t>> rd_paths_trunc(row_ids.size());
 
-    const graph::boss::BOSS &boss = graph_->get_boss();
-    const bit_vector &rd_succ = fork_succ_.size() ? fork_succ_ : boss.get_last();
-
     for (size_t i = 0; i < row_ids.size(); ++i) {
         Row row = row_ids[i];
 
-        graph::boss::BOSS::edge_index boss_edge = graph_->kmer_to_boss_index(
-                graph::AnnotatedSequenceGraph::anno_to_graph_index(row));
-
         while (true) {
-            row = graph::AnnotatedSequenceGraph::graph_to_anno_index(
-                    graph_->boss_to_kmer_index(boss_edge));
-
             auto [it, is_new] = node_to_rd.try_emplace(row, rd_ids.size());
             rd_paths_trunc[i].push_back(it.value());
 
@@ -230,7 +215,9 @@ RowDiff<BaseMatrix>::get_rows(const std::vector<Row> &row_ids) const {
             if (anchor_[row])
                 break;
 
-            boss_edge = boss.row_diff_successor(boss_edge, rd_succ);
+            auto node = graph::AnnotatedSequenceGraph::anno_to_graph_index(row);
+            node = graph_->row_diff_successor(node, fork_succ_);
+            row = graph::AnnotatedSequenceGraph::graph_to_anno_index(node);
         }
     }
 

metagraph/src/annotation/int_matrix/row_diff/int_row_diff.hpp

@@ -13,7 +13,6 @@
 #include "common/logger.hpp"
 #include "common/utils/template_utils.hpp"
 #include "graph/annotated_dbg.hpp"
-#include "graph/representation/succinct/boss.hpp"
 #include "graph/representation/succinct/dbg_succinct.hpp"
 #include "annotation/binary_matrix/row_diff/row_diff.hpp"
 #include "annotation/int_matrix/base/int_matrix.hpp"
@@ -86,9 +85,9 @@ template <class BaseMatrix>
 std::vector<IntMatrix::Row> IntRowDiff<BaseMatrix>::get_column(Column j) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     const graph::boss::BOSS &boss = graph_->get_boss();
-    assert(!fork_succ_.size() || fork_succ_.size() == boss.get_last().size());
 
     // TODO: implement a more efficient algorithm
     std::vector<Row> result;
@@ -113,14 +112,15 @@ std::vector<IntMatrix::RowValues>
 IntRowDiff<BaseMatrix>::get_row_values(const std::vector<Row> &row_ids) const {
     assert(graph_ && "graph must be loaded");
     assert(anchor_.size() == diffs_.num_rows() && "anchors must be loaded");
-    assert(!fork_succ_.size() || fork_succ_.size() == graph_->get_boss().get_last().size());
+    assert(!fork_succ_.size() || fork_succ_.size() == graph_->num_nodes() + 1);
 
     // get row-diff paths
     auto [rd_ids, rd_paths_trunc] = get_rd_ids(row_ids);
 
     std::vector<RowValues> rd_rows = diffs_.get_row_values(rd_ids);
     for (auto &row : rd_rows) {
         decode_diffs(&row);
+        std::sort(row.begin(), row.end());
     }
 
     rd_ids = std::vector<Row>();
@@ -129,17 +129,15 @@ IntRowDiff<BaseMatrix>::get_row_values(const std::vector<Row> &row_ids) const {
     std::vector<RowValues> rows(row_ids.size());
 
     for (size_t i = 0; i < row_ids.size(); ++i) {
-        RowValues &result = rows[i];
+        const auto &rd_path = rd_paths_trunc[i];
         // propagate back and reconstruct full annotations for predecessors
-        for (auto it = rd_paths_trunc[i].rbegin(); it != rd_paths_trunc[i].rend(); ++it) {
-            std::sort(rd_rows[*it].begin(), rd_rows[*it].end());
-            add_diff(rd_rows[*it], &result);
-            // replace diff row with full reconstructed annotation
-            rd_rows[*it] = result;
+        for (size_t j = 0; j + 1 < rd_path.size(); ++j) {
+            auto [node, succ] = rd_path[j];
+            // reconstruct annotation by adding the diff (full succ + diff)
+            add_diff(rd_rows[succ], &rd_rows[node]);
         }
-        assert(std::all_of(result.begin(), result.end(),
-                           [](auto &p) { return p.second; }));
-        assert(std::all_of(result.begin(), result.end(),
+        rows[i] = rd_rows[rd_path.back().second];
+        assert(std::all_of(rows[i].begin(), rows[i].end(),
                            [](auto &p) { return (int64_t)p.second > 0; }));
     }