graph.hpp

Go to the documentation of this file.

#pragma once

 
#include "implementation/lattice.hpp"
 
#include <algorithm>
#include <array>
#include <cassert>
#include <concepts>
#include <cstddef>
#include <iterator>
#include <optional>
#include <queue>
#include <ranges>
#include <stack>
#include <type_traits>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
 
namespace pgl {

template <class Vertex>
class Graph {
    using AdjacencyMap = std::unordered_map<Vertex, std::unordered_set<Vertex>>;
 
public:
    using VertexType = Vertex;

    using NeighborSet = std::unordered_set<Vertex>;

    using EdgeType = std::array<Vertex, 2>;

    class Iterator {
        using BaseIterator = typename AdjacencyMap::const_iterator;
 
    public:
        using iterator_category = std::forward_iterator_tag;
        using iterator_concept = std::forward_iterator_tag;
        using difference_type = std::ptrdiff_t;
        using value_type = Vertex;
        using pointer = const Vertex*;
        using reference = const Vertex&;

        Iterator() = default;

        reference operator*() const {
            return iterator_->first;
        }

        pointer operator->() const {
            return &iterator_->first;
        }

        Iterator& operator++() {
            ++iterator_;
            return *this;
        }

        Iterator operator++(int) {
            Iterator previous = *this;
            ++(*this);
            return previous;
        }
 
        friend bool operator==(const Iterator&, const Iterator&) = default;
 
    private:
        friend class Graph;
 
        explicit Iterator(BaseIterator iterator)
            : iterator_(iterator) {}
 
        BaseIterator iterator_;
    };

    class EdgeIterator {
        using OuterIterator = typename AdjacencyMap::const_iterator;
        using InnerIterator = typename NeighborSet::const_iterator;
 
    public:
        // Multi-pass, hence a forward iterator to the C++20 concepts, but the
        // edge is built on dereference rather than stored, so the old category
        // stays "input": only the concept admits a reference that is a value.
        using iterator_category = std::input_iterator_tag;
        using iterator_concept = std::forward_iterator_tag;
        using difference_type = std::ptrdiff_t;
        using value_type = EdgeType;
        using reference = EdgeType;

        EdgeIterator() = default;

        reference operator*() const {
            return EdgeType{outer_->first, *inner_};
        }

        EdgeIterator& operator++() {
            ++inner_;
            skipToEdge();
            return *this;
        }

        EdgeIterator operator++(int) {
            EdgeIterator previous = *this;
            ++(*this);
            return previous;
        }
 
        friend bool operator==(const EdgeIterator&, const EdgeIterator&) = default;
 
    private:
        friend class Graph;
 
        EdgeIterator(OuterIterator iterator, OuterIterator last)
            : outer_(iterator), last_(last) {
            if (outer_ != last_) {
                inner_ = outer_->second.cbegin();
                skipToEdge();
            }
        }
 
        // Advances the position to the next neighbor larger than its own
        // vertex, which is the one direction of an edge that this iterator
        // reports. Past the last one, the inner iterator is reset to its
        // value-initialized state so that it compares equal to the end
        // iterator's, whatever adjacency set the walk stopped in.
        void skipToEdge() {
            while (outer_ != last_) {
                if (inner_ == outer_->second.cend()) {
                    ++outer_;
                    if (outer_ == last_) {
                        break;
                    }
                    inner_ = outer_->second.cbegin();
                } else if (outer_->first < *inner_) {
                    return;
                } else {
                    ++inner_;
                }
            }
            inner_ = InnerIterator{};
        }
 
        OuterIterator outer_;
        OuterIterator last_;
        InnerIterator inner_;
    };
 
    using iterator = Iterator;
    using const_iterator = Iterator;

    Graph() = default;

    explicit Graph(const std::vector<std::array<Vertex, 2>>& edges) {
        for (const auto& [u, v] : edges) {
            addEdge(u, v);
        }
    }

    void addVertex(const Vertex& vertex) {
        adjacency_.try_emplace(vertex);
    }

    void addEdge(const Vertex& u, const Vertex& v) {
        if (u != v) {
            adjacency_[u].insert(v);
            adjacency_[v].insert(u);
        }
    }

    [[nodiscard]] bool containsVertex(const Vertex& vertex) const {
        return adjacency_.contains(vertex);
    }

    [[nodiscard]] bool containsEdge(const Vertex& u, const Vertex& v) const {
        const auto uIt = adjacency_.find(u);
        return uIt != adjacency_.end() && uIt->second.contains(v);
    }

    [[nodiscard]] int degree(const Vertex& vertex) const {
        const auto vertexIt = adjacency_.find(vertex);
        if (vertexIt == adjacency_.end()) {
            return -1;
        }
        return static_cast<int>(vertexIt->second.size());
    }

    [[nodiscard]] int maxDegree() const {
        int result = -1;
        for (const auto& entry : adjacency_) {
            result = std::max(result, static_cast<int>(entry.second.size()));
        }
        return result;
    }

    [[nodiscard]] int vertexCount() const {
        return static_cast<int>(adjacency_.size());
    }

    [[nodiscard]] int edgeCount() const {
        std::size_t directedEdgeCount = 0;
        for (const auto& entry : adjacency_) {
            directedEdgeCount += entry.second.size();
        }
        assert(directedEdgeCount % 2 == 0);
        return static_cast<int>(directedEdgeCount / 2);
    }

    void removeEdge(const Vertex& u, const Vertex& v) {
        const auto uIt = adjacency_.find(u);
        if (uIt == adjacency_.end() || !uIt->second.contains(v)) {
            return;
        }
        uIt->second.erase(v);
        adjacency_.at(v).erase(u);
    }

    void removeVertex(const Vertex& vertex) {
        const auto vertexIt = adjacency_.find(vertex);
        if (vertexIt == adjacency_.end()) {
            return;
        }
 
        // Copy because erasing incident edges invalidates adjacency iterators.
        const NeighborSet neighbors = vertexIt->second;
        for (const Vertex& neighbor : neighbors) {
            adjacency_.at(neighbor).erase(vertex);
        }
        adjacency_.erase(vertexIt);
    }

    void clear() {
        adjacency_.clear();
    }

    [[nodiscard]] auto vertices() const {
        return std::ranges::subrange(begin(), end());
    }

    [[nodiscard]] auto edges() const
        requires std::totally_ordered<Vertex>
    {
        return std::ranges::subrange(
            EdgeIterator(adjacency_.cbegin(), adjacency_.cend()),
            EdgeIterator(adjacency_.cend(), adjacency_.cend())
        );
    }

    [[nodiscard]] const NeighborSet& neighbors(const Vertex& vertex) const {
        return adjacency_.at(vertex);
    }

    [[nodiscard]] NeighborSet closedNeighbors(const Vertex& vertex) const {
        NeighborSet result = adjacency_.at(vertex);
        result.insert(vertex);
        return result;
    }

    [[nodiscard]] std::vector<Vertex> bfs(const Vertex& vertex, int maxVertices = 0) const {
        if (!containsVertex(vertex) || maxVertices < 0) {
            return {};
        }
 
        const std::size_t limit = maxVertices == 0
            ? adjacency_.size()
            : static_cast<std::size_t>(maxVertices);
 
        NeighborSet visited;
        std::vector<Vertex> result;
        std::queue<Vertex> queue;
 
        visited.insert(vertex);
        queue.push(vertex);
        while (!queue.empty() && result.size() < limit) {
            const Vertex current = queue.front();
            queue.pop();
            result.push_back(current);
 
            for (const Vertex& neighbor : neighbors(current)) {
                if (visited.insert(neighbor).second) {
                    queue.push(neighbor);
                }
            }
        }
 
        return result;
    }

    [[nodiscard]] std::vector<std::vector<Vertex>> components() const {
        NeighborSet visited;
        std::vector<std::vector<Vertex>> result;
 
        for (const auto& entry : adjacency_) {
            if (visited.contains(entry.first)) {
                continue;
            }
 
            result.push_back(bfs(entry.first));
            visited.insert(result.back().begin(), result.back().end());
        }
 
        std::sort(result.begin(), result.end(), [](const auto& a, const auto& b) {
            return a.size() > b.size();
        });
        return result;
    }

    [[nodiscard]] std::vector<std::vector<Vertex>> biconnectedComponents() const {
        struct SearchData {
            std::size_t index;
            std::size_t low;
        };
 
        struct Frame {
            Vertex vertex;
            std::optional<Vertex> parent;
            typename NeighborSet::const_iterator nextNeighbor;
            typename NeighborSet::const_iterator endNeighbor;
        };
 
        std::size_t timer = 0;
        std::unordered_map<Vertex, SearchData> searchData;
        std::vector<std::pair<Vertex, Vertex>> edgeStack;
        std::vector<std::vector<Vertex>> result;
 
        for (const auto& entry : adjacency_) {
            const Vertex& start = entry.first;
            if (searchData.contains(start)) {
                continue;
            }
 
            searchData.emplace(start, SearchData{timer, timer});
            ++timer;
 
            std::stack<Frame> stack;
            stack.push(Frame{start, std::nullopt, entry.second.cbegin(), entry.second.cend()});
 
            while (!stack.empty()) {
                Frame& frame = stack.top();
                const Vertex vertex = frame.vertex;
 
                if (frame.nextNeighbor != frame.endNeighbor) {
                    const Vertex neighbor = *frame.nextNeighbor;
                    ++frame.nextNeighbor;
 
                    if (frame.parent.has_value() && neighbor == *frame.parent) {
                        continue;
                    }
 
                    const auto neighborData = searchData.find(neighbor);
                    if (neighborData == searchData.end()) {
                        edgeStack.emplace_back(vertex, neighbor);
                        searchData.emplace(neighbor, SearchData{timer, timer});
                        ++timer;
 
                        const NeighborSet& nextNeighbors = adjacency_.at(neighbor);
                        stack.push(Frame{
                            neighbor,
                            vertex,
                            nextNeighbors.cbegin(),
                            nextNeighbors.cend(),
                        });
                    } else if (neighborData->second.index < searchData.at(vertex).index) {
                        SearchData& vertexData = searchData.at(vertex);
                        vertexData.low = std::min(vertexData.low, neighborData->second.index);
                        edgeStack.emplace_back(vertex, neighbor);
                    }
                    continue;
                }
 
                const std::optional<Vertex> parent = frame.parent;
                stack.pop();
                if (!parent.has_value()) {
                    continue;
                }
 
                SearchData& parentData = searchData.at(*parent);
                const SearchData& vertexData = searchData.at(vertex);
                parentData.low = std::min(parentData.low, vertexData.low);
 
                if (vertexData.low >= parentData.index) {
                    NeighborSet addedVertices;
                    std::vector<Vertex> component;
                    bool foundTreeEdge = false;
 
                    do {
                        assert(!edgeStack.empty());
                        const auto edge = std::move(edgeStack.back());
                        edgeStack.pop_back();
 
                        if (addedVertices.insert(edge.first).second) {
                            component.push_back(edge.first);
                        }
                        if (addedVertices.insert(edge.second).second) {
                            component.push_back(edge.second);
                        }
                        foundTreeEdge = edge.first == *parent && edge.second == vertex;
                    } while (!foundTreeEdge);
 
                    result.push_back(std::move(component));
                }
            }
        }
 
        std::sort(result.begin(), result.end(), [](const auto& a, const auto& b) {
            return a.size() > b.size();
        });
        return result;
    }

    [[nodiscard]] std::vector<std::vector<Vertex>> cliqueCover() const {
        // DSATUR works with vertex indices, so the vertices are materialized
        // once rather than taken from the lazy view.
        const std::vector<Vertex> graphVertices(begin(), end());
        const std::size_t vertexCount = graphVertices.size();
        const std::size_t uncolored = vertexCount;
 
        std::vector<std::size_t> colors(vertexCount, uncolored);
        std::vector<std::size_t> complementDegrees(vertexCount);
        std::vector<std::unordered_set<std::size_t>> neighborColors(vertexCount);
        std::vector<std::vector<Vertex>> result;
 
        for (std::size_t i = 0; i < vertexCount; ++i) {
            complementDegrees[i] = vertexCount - 1 - adjacency_.at(graphVertices[i]).size();
        }
 
        for (std::size_t coloredCount = 0; coloredCount < vertexCount; ++coloredCount) {
            std::size_t selected = uncolored;
            for (std::size_t i = 0; i < vertexCount; ++i) {
                if (colors[i] != uncolored) {
                    continue;
                }
 
                if (selected == uncolored ||
                    neighborColors[i].size() > neighborColors[selected].size() ||
                    (neighborColors[i].size() == neighborColors[selected].size() &&
                     complementDegrees[i] > complementDegrees[selected])) {
                    selected = i;
                }
            }
            assert(selected != uncolored);
 
            std::size_t color = 0;
            while (neighborColors[selected].contains(color)) {
                ++color;
            }
            colors[selected] = color;
 
            if (color == result.size()) {
                result.emplace_back();
            }
            result[color].push_back(graphVertices[selected]);
 
            const NeighborSet& selectedNeighbors = adjacency_.at(graphVertices[selected]);
            for (std::size_t i = 0; i < vertexCount; ++i) {
                if (colors[i] == uncolored && i != selected &&
                    !selectedNeighbors.contains(graphVertices[i])) {
                    neighborColors[i].insert(color);
                }
            }
        }
 
        std::sort(result.begin(), result.end(), [](const auto& a, const auto& b) {
            return a.size() > b.size();
        });
        return result;
    }

    [[nodiscard]] std::vector<Vertex> independentSet() const {
        std::vector<Vertex> orderedVertices(begin(), end());
        std::sort(orderedVertices.begin(), orderedVertices.end(), [this](
            const Vertex& a,
            const Vertex& b
        ) {
            return adjacency_.at(a).size() < adjacency_.at(b).size();
        });
 
        NeighborSet selected;
        std::vector<Vertex> result;
        for (const Vertex& vertex : orderedVertices) {
            const NeighborSet& vertexNeighbors = adjacency_.at(vertex);
            if (std::ranges::none_of(vertexNeighbors, [&selected](const Vertex& neighbor) {
                    return selected.contains(neighbor);
                })) {
                selected.insert(vertex);
                result.push_back(vertex);
            }
        }
        return result;
    }

    template <class WeightFunction>
    [[nodiscard]] Graph spanningTree(WeightFunction weight) const {
        using Weight = std::invoke_result_t<WeightFunction&, const Vertex&, const Vertex&>;
 
        struct Candidate {
            Weight weight;
            Vertex from;
            Vertex to;
        };
 
        // A priority_queue pops its largest element, so order edges by
        // decreasing weight to obtain the lightest frontier edge.
        const auto heavier = [](const Candidate& a, const Candidate& b) {
            return b.weight < a.weight;
        };
 
        Graph result;
        NeighborSet visited;
        std::priority_queue<Candidate, std::vector<Candidate>, decltype(heavier)> frontier(heavier);
 
        const auto pushIncidentEdges = [&](const Vertex& vertex) {
            for (const Vertex& neighbor : adjacency_.at(vertex)) {
                if (!visited.contains(neighbor)) {
                    frontier.push(Candidate{weight(vertex, neighbor), vertex, neighbor});
                }
            }
        };
 
        for (const auto& entry : adjacency_) {
            if (visited.contains(entry.first)) {
                continue;
            }
 
            // Grow one tree per connected component. The frontier is empty
            // again once a component is exhausted, so it can be reused.
            visited.insert(entry.first);
            result.addVertex(entry.first);
            pushIncidentEdges(entry.first);
 
            while (!frontier.empty()) {
                const Candidate best = frontier.top();
                frontier.pop();
                if (!visited.insert(best.to).second) {
                    continue;
                }
                result.addEdge(best.from, best.to);
                pushIncidentEdges(best.to);
            }
        }
 
        return result;
    }

    template <class WeightFunction>
    [[nodiscard]] std::vector<Vertex> shortestPath(
        const Vertex& source,
        const Vertex& target,
        WeightFunction weight
    ) const {
        using Weight = std::invoke_result_t<WeightFunction&, const Vertex&, const Vertex&>;
 
        struct Candidate {
            Weight weight;
            Vertex from;
            Vertex to;
        };
 
        if (!containsVertex(source) || !containsVertex(target)) {
            return {};
        }
        if (source == target) {
            return {source};
        }
 
        // A priority_queue pops its largest element, so order candidates by
        // decreasing distance to obtain the closest unsettled vertex.
        const auto farther = [](const Candidate& a, const Candidate& b) {
            return b.weight < a.weight;
        };
 
        NeighborSet settled;
        std::unordered_map<Vertex, Vertex> parent;
        std::priority_queue<Candidate, std::vector<Candidate>, decltype(farther)> frontier(farther);
 
        const auto pushIncidentEdges = [&](const Vertex& vertex, const Candidate& candidate) {
            for (const Vertex& neighbor : adjacency_.at(vertex)) {
                if (!settled.contains(neighbor)) {
                    frontier.push(Candidate{
                        candidate.weight + weight(vertex, neighbor),
                        vertex,
                        neighbor,
                    });
                }
            }
        };
 
        // Seeding the frontier with the edges leaving the source, rather than
        // with the source itself, keeps every distance a sum of edge weights:
        // the weight type needs no zero of its own.
        settled.insert(source);
        for (const Vertex& neighbor : adjacency_.at(source)) {
            frontier.push(Candidate{weight(source, neighbor), source, neighbor});
        }
 
        while (!frontier.empty()) {
            const Candidate best = frontier.top();
            frontier.pop();
            if (!settled.insert(best.to).second) {
                continue;
            }
            parent.emplace(best.to, best.from);
 
            if (best.to == target) {
                std::vector<Vertex> result{target};
                while (result.back() != source) {
                    result.push_back(parent.at(result.back()));
                }
                std::reverse(result.begin(), result.end());
                return result;
            }
 
            pushIncidentEdges(best.to, best);
        }
 
        return {};
    }

    template <class WeightFunction, class LowerBoundFunction>
    [[nodiscard]] std::vector<Vertex> shortestPath(
        const Vertex& source,
        const Vertex& target,
        WeightFunction weight,
        LowerBoundFunction lowerBound
    ) const {
        using Weight =
            std::invoke_result_t<WeightFunction&, const Vertex&, const Vertex&>;
 
        struct Candidate {
            Weight distance;
            Weight estimate;
            Vertex vertex;
        };
 
        if (!containsVertex(source) || !containsVertex(target)) {
            return {};
        }
        if (source == target) {
            return {source};
        }
 
        // A priority_queue pops its largest element, so order candidates by
        // decreasing estimated total distance to obtain the most promising
        // unsettled vertex.
        const auto farther = [](const Candidate& a, const Candidate& b) {
            return b.estimate < a.estimate;
        };
 
        std::unordered_map<Vertex, Weight> distance;
        std::unordered_map<Vertex, Vertex> parent;
        std::priority_queue<Candidate, std::vector<Candidate>, decltype(farther)>
            frontier(farther);
 
        const auto relax = [&](const Vertex& from,
                               const Vertex& to,
                               const Weight& fromDistance) {
            if (to == source) {
                return;
            }
 
            const Weight newDistance = fromDistance + weight(from, to);
            const auto known = distance.find(to);
            if (known != distance.end() && !(newDistance < known->second)) {
                return;
            }
 
            distance.insert_or_assign(to, newDistance);
            parent.insert_or_assign(to, from);
            frontier.push(Candidate{
                newDistance,
                newDistance + lowerBound(to, target),
                to,
            });
        };
 
        // As in the Dijkstra overload, seed with the source's edges so the
        // weight type does not need a default-constructed zero.
        for (const Vertex& neighbor : adjacency_.at(source)) {
            const Weight neighborDistance = weight(source, neighbor);
            distance.emplace(neighbor, neighborDistance);
            parent.emplace(neighbor, source);
            frontier.push(Candidate{
                neighborDistance,
                neighborDistance + lowerBound(neighbor, target),
                neighbor,
            });
        }
 
        while (!frontier.empty()) {
            const Candidate best = frontier.top();
            frontier.pop();
 
            const auto known = distance.find(best.vertex);
            if (known == distance.end() || known->second < best.distance) {
                continue;
            }
 
            if (best.vertex == target) {
                std::vector<Vertex> result{target};
                while (result.back() != source) {
                    result.push_back(parent.at(result.back()));
                }
                std::reverse(result.begin(), result.end());
                return result;
            }
 
            for (const Vertex& neighbor : adjacency_.at(best.vertex)) {
                relax(best.vertex, neighbor, best.distance);
            }
        }
 
        return {};
    }

    [[nodiscard]] iterator begin() {
        return iterator(adjacency_.cbegin());
    }

    [[nodiscard]] iterator end() {
        return iterator(adjacency_.cend());
    }

    [[nodiscard]] const_iterator begin() const {
        return const_iterator(adjacency_.cbegin());
    }

    [[nodiscard]] const_iterator end() const {
        return const_iterator(adjacency_.cend());
    }

    [[nodiscard]] const_iterator cbegin() const {
        return begin();
    }

    [[nodiscard]] const_iterator cend() const {
        return end();
    }
 
private:
    AdjacencyMap adjacency_;
};
 
}  // namespace pgl