BlossomHelper.h

Go to the documentation of this file.

 
#pragma once
 
#include <ogdf/basic/EpsilonTest.h>
#include <ogdf/basic/Graph.h>
#include <ogdf/basic/GraphCopy.h>
#include <ogdf/basic/GraphList.h>
#include <ogdf/basic/Logger.h>
#include <ogdf/basic/basic.h>
#include <ogdf/graphalg/matching_blossom/Cycle.h>
#include <ogdf/graphalg/matching_blossom/Pseudonode.h>
#include <ogdf/graphalg/matching_blossom/utils.h>
 
#include <array>
#include <cstddef>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <unordered_map>
#include <unordered_set>
#include <utility>
#include <vector>
 
namespace ogdf {
namespace matching_blossom {
 
template<class TWeight>
class BlossomHelper : private Logger {
protected:
    bool m_greedyInit;
 
    GraphCopySimple m_graph;
 
    EdgeArray<TWeight> m_c;
    NodeArray<TWeight> m_y;
 
    NodeArray<edge> m_matching;
 
    std::unordered_map<node, Pseudonode*> m_pseudonodes;
 
    std::vector<node> m_repr;
 
    std::vector<node> m_shortcuts;
 
    EpsilonTest m_eps;
 
    // Multiply all integer weights by 4 to ensure that all primal and dual values stay integral.
    static const int WEIGHT_FACTOR = std::numeric_limits<TWeight>::is_integer ? 4 : 1;
 
    bool initDualSolution(NodeArray<TWeight>& minY) {
        // for _some_ reason, the array has to be initialized with 0, otherwise the algorithm fails
        m_y.init(m_graph, 0);
        for (node v : m_graph.nodes) {
            if (v->degree() == 0) {
                return false;
            }
            m_y[v] = minY[v];
        }
        if (!m_greedyInit) {
            return true;
        }
        for (node v : m_graph.nodes) {
            if (matching(v) != nullptr) {
                continue;
            }
            edge matchingEdge = nullptr;
            TWeight maxChange = infinity<TWeight>();
            for (auto adj : v->adjEntries) {
                edge e = adj->theEdge();
                node w = adj->twinNode();
                TWeight reducedWeight = getReducedWeight(e);
                if (m_eps.leq(reducedWeight, maxChange)) {
                    maxChange = reducedWeight;
                    if (!matching(w)) {
                        matchingEdge = e;
                    }
                }
            }
            y(v) += maxChange;
            if (matchingEdge != nullptr && isEqualityEdge(matchingEdge)) {
                addToMatching(matchingEdge);
            }
        }
#ifdef OGDF_DEBUG
        for (node v : m_graph.nodes) {
            lout() << "initial y for node " << v << ": " << y(v) << std::endl;
            edge e = matching(v);
            if (e && e->source() == v) {
                OGDF_ASSERT(e == matching(e->target()));
                OGDF_ASSERT(isEqualityEdge(e));
                lout() << "initial matching edge: " << e << std::endl;
            }
        }
        for (edge e : m_graph.edges) {
            OGDF_ASSERT(getReducedWeight(e) >= 0);
            if (isEqualityEdge(e)) {
                OGDF_ASSERT(matching(e->source()) || matching(e->target()));
            }
        }
#endif
        return true;
    }
 
    node findParentInRepr(node v, node child = nullptr) {
        node parent = m_repr[v->index()];
        OGDF_ASSERT(parent != nullptr);
        if (parent == v) {
            return child;
        } else {
            node grandparent = findParentInRepr(parent, v);
            m_shortcuts[v->index()] = grandparent;
            return grandparent;
        }
    }
 
    void deletePseudonodes() {
        for (auto entry : m_pseudonodes) {
            delete entry.second;
        }
        m_pseudonodes.clear();
    }
 
public:
#ifdef OGDF_HEAVY_DEBUG
    void getReducedEdgeWeights(std::unordered_map<edge, TWeight>& edges) {
        for (edge e : m_graph.edges) {
            if (repr(e->source()) == e->source() && repr(e->target()) == e->target()) {
                edges[e] = getRealReducedWeight(e);
            }
        }
    }
 
    void assertConsistentEdgeWeights(std::unordered_map<edge, TWeight>& edges) {
        for (edge e : m_graph.edges) {
            if (repr(e->source()) == e->source() && repr(e->target()) == e->target()
                    && edges.find(e) != edges.end()) {
                OGDF_ASSERT(edges[e] == getRealReducedWeight(e));
            }
        }
    }
#endif
 
    BlossomHelper(bool greedyInit) : m_greedyInit(greedyInit) { }
 
    ~BlossomHelper() { deletePseudonodes(); }
 
    template<class WeightContainer>
    bool init(const Graph& graph, const WeightContainer& weights) {
        if (graph.numberOfNodes() % 2 == 1) {
            return false;
        }
        m_graph.clear();
        m_graph.init(graph);
        m_c.init(m_graph);
        NodeArray<TWeight> minY(m_graph, -1);
        for (edge e : m_graph.edges) {
            if (e->isSelfLoop()) {
                throw std::runtime_error("Self-loops are not supported.");
            }
            // copy initial edge costs from original edge, multiplied by WEIGHT_FACTOR
            TWeight weight = WEIGHT_FACTOR * getWeight<TWeight>(m_graph.original(e), weights);
            m_c[e] = weight;
            // init dual variables for all nodes
            auto halfWeight = weight / 2;
            for (node v : e->nodes()) {
                if (minY[v] == -1 || m_eps.less(halfWeight, minY[v])) {
                    minY[v] = halfWeight;
                }
            }
        }
        m_matching.init(m_graph, nullptr);
        deletePseudonodes();
        m_repr.clear();
        m_shortcuts.clear();
        for (node v : m_graph.nodes) {
            m_repr.push_back(v);
            m_shortcuts.push_back(v);
        }
        return initDualSolution(minY);
    }
 
    /* Getters */
 
    GraphCopySimple& graph() { return m_graph; }
 
    TWeight& c(edge e) { return m_c[e]; }
 
    TWeight& y(node v) { return m_y[v]; }
 
    NodeArray<edge>& matching() { return m_matching; }
 
    edge& matching(node v) { return m_matching[v]; }
 
    std::unordered_map<node, Pseudonode*>& pseudonodes() { return m_pseudonodes; }
 
    Pseudonode* pseudonode(node v) { return tryGetPointerFromMap(m_pseudonodes, v); }
 
    bool isPseudonode(node v) { return pseudonode(v) != nullptr; }
 
    /* End of getters */
 
    node repr(node v) {
        node parent = m_repr[v->index()];
        if (parent == v || parent == nullptr) {
            return v;
        }
        node child = reprChild(v);
        parent = m_repr[child->index()];
        if (parent == nullptr) {
            return child;
        } else {
            return parent;
        }
    }
 
    node reprChild(node v) {
        node shortcut = m_shortcuts[v->index()];
        node parent = m_repr[shortcut->index()];
        if (parent == shortcut || parent == nullptr) {
            return findParentInRepr(v);
        } else {
            parent = findParentInRepr(shortcut);
            m_shortcuts[v->index()] = parent;
            return parent;
        }
    }
 
    void addPseudonode(Pseudonode* pseudonode) {
        m_pseudonodes[pseudonode->graphNode] = pseudonode;
        m_repr.push_back(pseudonode->graphNode);
        m_shortcuts.push_back(pseudonode->graphNode);
        OGDF_ASSERT(m_repr[pseudonode->graphNode->index()] == pseudonode->graphNode);
        for (node v : pseudonode->cycle->nodes()) {
            m_repr[v->index()] = pseudonode->graphNode;
        }
    }
 
    void removePseudonode(Pseudonode* pseudonode) {
        node theNode = pseudonode->graphNode;
        m_pseudonodes.erase(theNode);
        m_matching[theNode] = nullptr;
        m_repr[theNode->index()] = nullptr;
        m_graph.delNode(theNode);
        for (node v : pseudonode->cycle->nodes()) {
            m_repr[v->index()] = v;
        }
    }
 
    void expandRepr(Pseudonode* pseudonode) {
        m_repr[pseudonode->graphNode->index()] = nullptr;
        for (node v : pseudonode->cycle->nodes()) {
            m_repr[v->index()] = v;
        }
    }
 
    node getBaseNode(edge e, node v) { return getBaseNodes(e, v).first; }
 
    node getOppositeBaseNode(edge e, node v) { return getBaseNodes(e, v).second; }
 
    std::pair<node, node> getBaseNodes(edge e, node v = nullptr) {
        edge eOrig = m_graph.original(e);
        node source = m_graph.copy(eOrig->source());
        node target = m_graph.copy(eOrig->target());
        if (v != nullptr && repr(target) == v) {
            return {target, source};
        } else {
            OGDF_ASSERT(v == nullptr || repr(source) == v);
            return {source, target};
        }
    }
 
    bool isZero(TWeight x) { return m_eps.equal(x, (TWeight)0); }
 
    TWeight getReducedWeight(edge e) { return c(e) - y(e->source()) - y(e->target()); }
 
    virtual TWeight getRealReducedWeight(edge e) { return getReducedWeight(e); }
 
    virtual bool isEqualityEdge(edge e) { return isZero(getRealReducedWeight(e)); }
 
    bool isZeroCostNode(node v) { return isZero(y(v)); }
 
    void getOriginalMatching(std::unordered_set<edge>& matching) {
        // Extend matching in the copy of the graph with the edges in the contracted cycles
        std::vector<Pseudonode*> stack;
        for (auto entry : m_pseudonodes) {
            if (repr(entry.first) == entry.first) {
                stack.push_back(entry.second);
            }
        }
 
        while (!stack.empty()) {
            auto pseudonode = stack.back();
            stack.pop_back();
            node graphNode = pseudonode->graphNode;
            auto cycle = pseudonode->cycle;
            // Find the edge entering the cycle
            edge matchingEdge = m_matching[graphNode];
            // Find the node in the cycle where the matching edge enters
            node matchedNodeInner = getBaseNode(matchingEdge, graphNode);
            node matchedNode = reprChild(matchedNodeInner);
            // Set the matching edge of the node in case it is a pseudonode and also needs to be
            // expanded later. This breaks the invariant that m_matching always contains exactly the two
            // endpoints pointing to the same edge, but since this is the last step of the algorithm, it
            // doesn't matter.
            m_matching[matchedNode] = matchingEdge;
            auto startIndex = cycle->indexOf(matchedNode);
            auto edges = cycle->edgeOrder();
            // Iterate the edges and add every second one from the perspective of the matchedNode to the
            // matching. Since startIndex is the index of the edge _before_ the matchedNode, we start at
            // index 2.
            for (size_t i = 2; i < edges.size(); i += 2) {
                edge e = edges[(startIndex + i) % edges.size()];
                addToMatching(e);
            }
            // Also expand all children
            for (node v : cycle->nodes()) {
                if (auto child = this->pseudonode(v)) {
                    stack.push_back(child);
                }
            }
            expandRepr(pseudonode);
        }
        // Add all respective edges of the original graph to the output parameter
        for (edge e : m_matching) {
            matching.insert(m_graph.original(e));
        }
    }
 
    void addToMatching(edge e) { m_matching[e->source()] = m_matching[e->target()] = e; }
};
 
}
}