Splitter part 1

src/common/CMakeLists.txt

@@ -14,6 +14,7 @@ add_subdirectory(pipeline)
 add_subdirectory(sequence)
 add_subdirectory(assembly_graph)
 add_subdirectory(alignment)
+add_subdirectory(auxiliary_graphs)
 add_subdirectory(modules/path_extend)
 add_subdirectory(modules)
 add_subdirectory(paired_info)
@@ -29,7 +30,7 @@ add_library(common_modules
             adt/concurrent_dsu.cpp)
 
 target_link_libraries(common_modules library
-                      assembly_graph input pipeline sequence
-                      coverage_model paired_info path_extend
-                      stages mph_index modules
+                      assembly_graph auxiliary_graphs input
+                      pipeline sequence coverage_model
+                      paired_info path_extend stages mph_index modules
                       utils configs alignment)

src/common/auxiliary_graphs/CMakeLists.txt

@@ -0,0 +1,18 @@
+############################################################################
+# Copyright (c) 2019 Saint Petersburg State University
+# All Rights Reserved
+# See file LICENSE for details.
+############################################################################
+
+project(auxiliary_graphs CXX)
+
+add_library(auxiliary_graphs STATIC
+            contracted_graph/contracted_graph.cpp
+            contracted_graph/contracted_graph_builder.cpp
+            contracted_graph/contracted_graph_helper.cpp
+            contracted_graph/graph_condensation.cpp
+            contracted_graph/contracted_statistics.cpp
+            scaffold_graph/scaffold_vertex.cpp
+)
+
+target_link_libraries(auxiliary_graphs assembly_graph)

src/common/auxiliary_graphs/contracted_graph/contracted_graph.cpp

@@ -0,0 +1,203 @@
+//***************************************************************************
+//* Copyright (c) 2019 Saint Petersburg State University
+//* All Rights Reserved
+//* See file LICENSE for details.
+//***************************************************************************
+
+#include "contracted_graph.hpp"
+
+namespace contracted_graph {
+
+void AdjacencyMap::InsertPair(const AdjacencyMap::VertexId &vertex, const AdjacencyMap::ScaffoldVertex &edge) {
+    data_[vertex].insert(edge);
+}
+AdjacencyMap::const_iterator AdjacencyMap::begin() const {
+    return data_.begin();
+}
+AdjacencyMap::const_iterator AdjacencyMap::end() const {
+    return data_.end();
+}
+void AdjacencyMap::RemovePair(const VertexId &vertex, const AdjacencyMap::ScaffoldVertex &edge) {
+    data_.at(vertex).erase(edge);
+    if (data_.at(vertex).empty()) {
+        data_.erase(vertex);
+    }
+}
+bool AdjacencyMap::Contains(const VertexId &vertex, const AdjacencyMap::ScaffoldVertex &edge) {
+    auto vertex_entry = data_.find(vertex);
+    if (vertex_entry == data_.end()) {
+        return false;
+    }
+    return vertex_entry->second.find(edge) != vertex_entry->second.end();
+}
+bool AdjacencyMap::empty() const {
+    return data_.empty();
+}
+size_t AdjacencyMap::size() const {
+    return data_.size();
+}
+
+void ContractedGraph::InsertVertex(const ContractedGraph::VertexId &vertex) {
+    if (vertices_.insert(vertex).second) {
+        AdjacencyMap empty;
+        incoming_[vertex] = empty;
+        outcoming_[vertex] = empty;
+    }
+}
+void ContractedGraph::InsertEdge(const ContractedGraph::VertexId &head, const ContractedGraph::VertexId &tail,
+                                 const ContractedGraph::ScaffoldVertex &edge) {
+    VERIFY_DEV(vertices_.find(head) != vertices_.end());
+    VERIFY_DEV(vertices_.find(tail) != vertices_.end());
+    outcoming_[head].InsertPair(tail, edge);
+    incoming_[tail].InsertPair(head, edge);
+}
+ContractedGraph::const_entry_iterator ContractedGraph::in_entry_begin(const ContractedGraph::VertexId &vertex) const {
+    return incoming_.at(vertex).begin();
+}
+ContractedGraph::const_entry_iterator ContractedGraph::in_entry_end(const ContractedGraph::VertexId &vertex) const {
+    return incoming_.at(vertex).end();
+}
+adt::iterator_range<ContractedGraph::const_entry_iterator> ContractedGraph::IncomingEntries(
+        const ContractedGraph::VertexId &vertex) const {
+    return adt::make_range(in_entry_begin(vertex), in_entry_end(vertex));
+}
+ContractedGraph::const_entry_iterator ContractedGraph::out_entry_begin(const ContractedGraph::VertexId &vertex) const {
+    return outcoming_.at(vertex).begin();
+}
+ContractedGraph::const_entry_iterator ContractedGraph::out_entry_end(const ContractedGraph::VertexId &vertex) const {
+    return outcoming_.at(vertex).end();
+}
+adt::iterator_range<ContractedGraph::const_entry_iterator> ContractedGraph::OutcomingEntries(
+        const ContractedGraph::VertexId &vertex) const {
+    return adt::make_range(out_entry_begin(vertex), out_entry_end(vertex));
+}
+size_t ContractedGraph::GetOutDegree(const ContractedGraph::VertexId &vertex) const {
+    size_t result = 0;
+    for (const auto &entry: outcoming_.at(vertex)) {
+        result += entry.second.size();
+    }
+    return result;
+}
+size_t ContractedGraph::GetInDegree(const ContractedGraph::VertexId &vertex) const {
+    size_t result = 0;
+    for (const auto &entry: incoming_.at(vertex)) {
+        result += entry.second.size();
+    }
+    return result;
+}
+size_t ContractedGraph::GetCapacity(const ContractedGraph::VertexId &vertex) const {
+    return capacity_.at(vertex);
+}
+void ContractedGraph::InsertCapacity(const ContractedGraph::VertexId &vertex, size_t capacity) {
+    capacity_[vertex] = capacity;
+}
+bool ContractedGraph::ContainsVertex(const ContractedGraph::VertexId &vertex) const {
+    return vertices_.find(vertex) != vertices_.end();
+}
+ContractedGraph::const_vertex_iterator ContractedGraph::begin() const {
+    return vertices_.begin();
+}
+ContractedGraph::const_vertex_iterator ContractedGraph::end() const {
+    return vertices_.end();
+}
+size_t ContractedGraph::size() const {
+    return vertices_.size();
+}
+size_t ContractedGraph::CountEdges() const {
+    size_t result = 0;
+    for (const auto &vertex: vertices()) {
+        result += GetOutDegree(vertex);
+    }
+    return result;
+}
+void ContractedGraph::RemoveEdge(const VertexId &head, const VertexId &tail, const ContractedGraph::ScaffoldVertex &edge) {
+    VERIFY_DEV(ContainsVertex(head));
+    VERIFY_DEV(ContainsVertex(tail));
+    auto &head_outcoming = outcoming_.at(head);
+    auto &tail_incoming = incoming_.at(tail);
+    if (not head_outcoming.Contains(tail, edge)) {
+        return;
+    }
+    VERIFY_DEV(tail_incoming.Contains(head, edge));
+    head_outcoming.RemovePair(tail, edge);
+    tail_incoming.RemovePair(head, edge);
+}
+ContractedGraph::ContractedGraph(const Graph &assembly_graph) : assembly_graph_(assembly_graph) {}
+
+const debruijn_graph::Graph &ContractedGraph::GetAssemblyGraph() const {
+    return assembly_graph_;
+}
+ContractedGraph::ScaffoldVertex ContractedGraph::conjugate(ContractedGraph::ScaffoldVertex edge) const {
+    return edge.GetConjugateFromGraph(assembly_graph_);
+}
+//std::string ContractedGraph::EdgeNucls(ContractedGraph::EdgeId edge) const {
+//    return edge.GetSequence(assembly_graph_);
+//}
+
+double ContractedGraph::coverage(ContractedGraph::EdgeId edge) const {
+    return edge.GetCoverageFromGraph(assembly_graph_);
+}
+size_t ContractedGraph::length(ContractedGraph::EdgeId edge) const {
+    return edge.GetLengthFromGraph(assembly_graph_);
+}
+size_t ContractedGraph::int_id(ContractedGraph::EdgeId edge) const {
+    return edge.int_id();
+}
+adt::iterator_range<ContractedGraph::const_vertex_iterator> ContractedGraph::vertices() const {
+    return adt::make_range(begin(), end());
+}
+ContractedGraph::const_edge_iterator ContractedGraph::in_edge_begin(const VertexId &vertex) const {
+    auto entry_begin = in_entry_begin(vertex);
+    if (not incoming_.at(vertex).empty()) {
+        return ContractedGraph::const_edge_iterator(entry_begin, entry_begin->second.begin(), in_entry_end(vertex));
+    }
+    return const_edge_iterator(entry_begin, empty_.end(), in_entry_end(vertex));
+}
+ContractedGraph::const_edge_iterator ContractedGraph::in_edge_end(const VertexId &vertex) const {
+    auto entry_end = in_entry_end(vertex);
+    auto entry_last = std::prev(entry_end);
+    if (not incoming_.at(vertex).empty()) {
+        return const_edge_iterator(entry_end, entry_last->second.end(), entry_end);
+    }
+    return const_edge_iterator(entry_end, empty_.end(), entry_end);
+}
+adt::iterator_range<ContractedGraph::const_edge_iterator> ContractedGraph::IncomingEdges(const VertexId &vertex) const {
+    return adt::make_range(in_edge_begin(vertex), in_edge_end(vertex));
+}
+
+ContractedGraph::const_edge_iterator ContractedGraph::out_edge_begin(const VertexId &vertex) const {
+    auto entry_begin = out_entry_begin(vertex);
+    if (not outcoming_.at(vertex).empty()) {
+        return ContractedGraph::const_edge_iterator(entry_begin, entry_begin->second.begin(), out_entry_end(vertex));
+    }
+    return const_edge_iterator(entry_begin, empty_.end(), out_entry_end(vertex));
+}
+ContractedGraph::const_edge_iterator ContractedGraph::out_edge_end(const VertexId &vertex) const {
+    auto entry_end = out_entry_end(vertex);
+    auto entry_last = std::prev(entry_end);
+    if (not outcoming_.at(vertex).empty()) {
+        return const_edge_iterator(entry_end, entry_last->second.end(), entry_end);
+    }
+    return const_edge_iterator(entry_end, empty_.end(), entry_end);
+}
+adt::iterator_range<ContractedGraph::const_edge_iterator> ContractedGraph::OutgoingEdges(const VertexId &vertex) const {
+    return adt::make_range(out_edge_begin(vertex), out_edge_end(vertex));
+}
+auto ContractedGraph::canonical_edges() const {
+    return assembly_graph_.canonical_edges();
+}
+ContractedGraph::VertexId ContractedGraph::conjugate(const ContractedGraph::VertexId &vertex) const {
+    return assembly_graph_.conjugate(vertex);
+}
+Sequence ContractedGraph::EdgeNucls(ContractedGraph::EdgeId edge) const {
+    VERIFY(edge.GetType() == scaffold_graph::ScaffoldVertexT::Edge);
+    assembly_graph_.EdgeNucls(edge.GetFirstEdge());
+}
+size_t ContractedGraph::IncomingEdgeCount(const ContractedGraph::VertexId &vertex) const {
+    return incoming_.at(vertex).size();
+}
+size_t ContractedGraph::OutgoingEdgeCount(const contracted_graph::ContractedGraph::VertexId &vertex) const {
+    return outcoming_.at(vertex).size();
+}
+
+}

src/common/auxiliary_graphs/contracted_graph/contracted_graph.hpp

@@ -0,0 +1,147 @@
+//***************************************************************************
+//* Copyright (c) 2019 Saint Petersburg State University
+//* All Rights Reserved
+//* See file LICENSE for details.
+//***************************************************************************
+
+#pragma once
+
+#include "auxiliary_graphs/scaffold_graph/scaffold_vertex.hpp"
+#include "adt/iterator_range.hpp"
+#include "assembly_graph/core/graph.hpp"
+
+namespace contracted_graph {
+class AdjacencyMap {
+ public:
+    typedef debruijn_graph::VertexId VertexId;
+    typedef debruijn_graph::EdgeId EdgeId;
+    typedef scaffold_graph::ScaffoldVertex ScaffoldVertex;
+    typedef std::map<VertexId, std::unordered_set<ScaffoldVertex>>::const_iterator const_iterator;
+    typedef std::map<VertexId, std::unordered_set<ScaffoldVertex>>::value_type value_type;
+
+    AdjacencyMap() = default;
+    AdjacencyMap(const VertexId &vertex, const ScaffoldVertex &edge) : data_({{vertex, {edge}}}) {}
+    void InsertPair(const VertexId &vertex, const ScaffoldVertex &edge);
+    void RemovePair(const VertexId &vertex, const ScaffoldVertex &edge);
+    bool Contains(const VertexId &vertex, const ScaffoldVertex &edge);
+    bool empty() const;
+    size_t size() const;
+
+    const_iterator begin() const;
+    const_iterator end() const;
+
+  private:
+    std::map<debruijn_graph::VertexId, std::unordered_set<ScaffoldVertex>> data_;
+};
+
+class ContractedGraph {
+ public:
+    typedef debruijn_graph::VertexId VertexId;
+    typedef debruijn_graph::Graph Graph;
+    typedef std::set<VertexId> VertexContainer;
+    typedef std::map<VertexId, AdjacencyMap> EdgeContainer;
+    typedef scaffold_graph::ScaffoldVertex ScaffoldVertex;
+    typedef AdjacencyMap::const_iterator const_entry_iterator;
+    typedef VertexContainer::const_iterator const_vertex_iterator;
+    typedef std::unordered_set<ScaffoldVertex>::const_iterator internal_edge_iterator;
+    typedef ScaffoldVertex EdgeId;
+
+    class const_edge_iterator : public boost::iterator_facade<const_edge_iterator,
+                                                              const ScaffoldVertex,
+                                                              boost::forward_traversal_tag> {
+      public:
+        explicit const_edge_iterator(const_entry_iterator entry_it,
+                                     internal_edge_iterator edge_it,
+                                     const_entry_iterator entry_end):
+            entry_it_(entry_it),
+            edge_it_(edge_it),
+            entry_end_(entry_end) {}
+
+      private:
+        friend class boost::iterator_core_access;
+
+        bool equal(const const_edge_iterator &other) const {
+            return entry_it_ == other.entry_it_ and edge_it_ == other.edge_it_ and entry_end_ == other.entry_end_;
+        }
+
+        const ScaffoldVertex &dereference() const {
+            return *edge_it_;
+        }
+
+        void increment() {
+            ++edge_it_;
+            if (edge_it_ == entry_it_->second.end()) {
+                ++entry_it_;
+                if (entry_it_ != entry_end_) {
+                    edge_it_ = entry_it_->second.begin();
+                }
+            }
+        }
+
+        const_entry_iterator entry_it_;
+        internal_edge_iterator edge_it_;
+        const_entry_iterator entry_end_;
+    };
+
+    explicit ContractedGraph(const Graph &assembly_graph);
+    virtual ~ContractedGraph() = default;
+    ContractedGraph(ContractedGraph &&other) = default;
+
+    void InsertVertex(const VertexId &vertex);
+    void InsertEdge(const VertexId &head, const VertexId &tail, const ScaffoldVertex &edge);
+    void RemoveEdge(const VertexId &head, const VertexId &tail, const ScaffoldVertex &edge);
+    size_t GetOutDegree(const VertexId &vertex) const;
+    size_t GetInDegree(const VertexId &vertex) const;
+    size_t GetCapacity(const VertexId &vertex) const;
+    void InsertCapacity(const VertexId &vertex, size_t capacity);
+    bool ContainsVertex(const VertexId &vertex) const;
+
+    const_entry_iterator in_entry_begin(const VertexId &vertex) const;
+    const_entry_iterator in_entry_end(const VertexId &vertex) const;
+    adt::iterator_range<const_entry_iterator> IncomingEntries(const VertexId &vertex) const;
+    const_entry_iterator out_entry_begin(const VertexId &vertex) const;
+    const_entry_iterator out_entry_end(const VertexId &vertex) const;
+    adt::iterator_range<const_entry_iterator> OutcomingEntries(const VertexId &vertex) const;
+
+    const_edge_iterator in_edge_begin(const VertexId &vertex) const;
+    const_edge_iterator in_edge_end(const VertexId &vertex) const;
+    adt::iterator_range<const_edge_iterator> IncomingEdges(const VertexId &vertex) const;
+    size_t IncomingEdgeCount(const VertexId &vertex) const;
+    const_edge_iterator out_edge_begin(const VertexId &vertex) const;
+    const_edge_iterator out_edge_end(const VertexId &vertex) const;
+    adt::iterator_range<const_edge_iterator> OutgoingEdges(const VertexId &vertex) const;
+    size_t OutgoingEdgeCount(const VertexId &vertex) const;
+
+    const_vertex_iterator begin() const;
+    const_vertex_iterator end() const;
+    adt::iterator_range<const_vertex_iterator> vertices() const;
+    size_t size() const;
+    size_t CountEdges() const;
+
+    //fixme also iterates over short edges
+    auto canonical_edges () const;
+
+    const Graph &GetAssemblyGraph() const;
+//    std::string EdgeNucls(EdgeId edge) const;
+//    std::string VertexNucls(VertexId vertex) const;
+    Sequence EdgeNucls(EdgeId edge) const;
+    double coverage(EdgeId edge) const;
+    size_t length(EdgeId edge) const;
+    size_t int_id(EdgeId edge) const;
+
+    ScaffoldVertex conjugate(ScaffoldVertex edge) const;
+    VertexId conjugate(const VertexId &vertex) const;
+
+ protected:
+    EdgeContainer outcoming_;
+    EdgeContainer incoming_;
+    VertexContainer vertices_;
+    std::map<VertexId, size_t> capacity_;
+
+    //for edge iterator
+    std::unordered_set<ScaffoldVertex> empty_;
+
+    //for compatibility with visualizers and stuff
+    const Graph &assembly_graph_;
+    };
+}