FullComponentGeneratorDreyfusWagner.h

Go to the documentation of this file.

 
#pragma once
 
#include <ogdf/basic/ArrayBuffer.h>
#include <ogdf/basic/Graph.h>
#include <ogdf/basic/Hashing.h>
#include <ogdf/basic/List.h>
#include <ogdf/basic/Math.h>
#include <ogdf/basic/SubsetEnumerator.h>
#include <ogdf/basic/basic.h>
 
#include <limits>
 
namespace ogdf {
template<typename T>
class EdgeWeightedGraph;
template<typename T>
class EdgeWeightedGraphCopy;
 
namespace steiner_tree {
 
template<typename T>
class FullComponentGeneratorDreyfusWagner {
    const EdgeWeightedGraph<T>& m_G; 
    const List<node>& m_terminals; 
    const NodeArray<bool>& m_isTerminal; 
    const NodeArray<NodeArray<T>>& m_distance; 
    const NodeArray<NodeArray<edge>>& m_pred; 
    SubsetEnumerator<node> m_terminalSubset; 
 
    using NodePairs = ArrayBuffer<NodePair>;
 
    struct DWMData {
        T cost;
        NodePairs nodepairs;
        ArrayBuffer<const DWMData*> subgraphs;
 
        DWMData(T _cost, NodePairs _nodepairs) : cost(_cost), nodepairs(_nodepairs) { }
 
        explicit DWMData(T _cost = std::numeric_limits<T>::max()) : cost(_cost) { }
 
        void invalidate() {
            nodepairs.clear();
            subgraphs.clear();
        }
 
        bool valid() const { return cost == 0 || !(nodepairs.empty() && subgraphs.empty()); }
 
        void add(const DWMData* other) {
            if (valid()) {
                if (other->valid()) {
                    subgraphs.push(other);
                } else {
                    invalidate();
                }
            }
            cost += other->cost;
        }
 
        void clear() {
            invalidate();
            cost = 0; // make it valid again
        }
 
        void add(NodePair nodepair, T c) {
            if (valid()) {
                nodepairs.push(nodepair);
            }
            cost += c;
        }
    };
 
    struct DWMSplit {
        T cost = std::numeric_limits<T>::max();
        const DWMData* subgraph1 = nullptr;
        const DWMData* subgraph2 = nullptr;
 
        void set(const DWMData* s1, const DWMData* s2) {
            subgraph1 = s1;
            subgraph2 = s2;
            cost = s1->cost + s2->cost;
        }
    };
 
    class SortedNodeListHashFunc;
 
    Hashing<List<node>, DWMData, SortedNodeListHashFunc> m_map;
 
    const DWMData* dataOf(const List<node>& key) const {
        OGDF_ASSERT(key.size() > 1);
        OGDF_ASSERT(m_map.member(key));
        return &m_map.lookup(key)->info();
    }
 
    T costOf(const List<node>& key) const {
        OGDF_ASSERT(key.size() > 1);
        if (key.size() == 2) { // a shortcut to avoid using the hash table
#ifdef OGDF_FULL_COMPONENT_GENERATION_TERMINAL_SSSP_AWARE
            if (m_isTerminal[key.front()]) {
                return m_distance[key.front()][key.back()];
            } else {
                OGDF_ASSERT(m_isTerminal[key.back()]);
                return m_distance[key.back()][key.front()];
            }
#else
            return m_distance[key.front()][key.back()];
#endif
        }
        return dataOf(key)->cost;
    }
 
    bool safeIfSumSmaller(const T summand1, const T summand2, const T compareValue) const {
#ifdef OGDF_FULL_COMPONENT_GENERATION_ALWAYS_SAFE
        return summand1 + summand2 < compareValue;
#else
        return summand1 < std::numeric_limits<T>::max() && summand2 < std::numeric_limits<T>::max()
                && summand1 + summand2 < compareValue;
#endif
    }
 
    static void sortedInserter(node w, List<node>& list, bool& inserted, node newNode) {
        if (!inserted && w->index() > newNode->index()) {
            list.pushBack(newNode);
            inserted = true;
        }
        list.pushBack(w);
    }
 
    void makeKey(List<node>& newSubset, node v) const {
        bool inserted = false;
        m_terminalSubset.forEachMember([&](node w) { sortedInserter(w, newSubset, inserted, v); });
        if (!inserted) {
            newSubset.pushBack(v);
        }
    }
 
    void makeKey(List<node>& newSubset, List<node>& newComplement,
            const SubsetEnumerator<node>& subset, node v) const {
        bool insertedIntoSubset = false;
        bool insertedIntoComplement = false;
        // Interestingly std::bind is much slower than using lambdas (at least on g++ 6.3)
        subset.forEachMemberAndNonmember(
                [&](node w) { sortedInserter(w, newSubset, insertedIntoSubset, v); },
                [&](node w) { sortedInserter(w, newComplement, insertedIntoComplement, v); });
        if (!insertedIntoSubset) {
            newSubset.pushBack(v);
        }
        if (!insertedIntoComplement) {
            newComplement.pushBack(v);
        }
    }
 
    void computeSplit(NodeArray<DWMSplit>& split, node v, SubsetEnumerator<node>& subset) const {
        if (split[v].subgraph1 != nullptr) { // already computed
            return;
        }
 
        DWMSplit& best = split[v];
        for (subset.begin(1, subset.numberOfMembersAndNonmembers() / 2); subset.valid();
                subset.next()) {
            List<node> newSubset, newComplement;
            makeKey(newSubset, newComplement, subset, v);
 
            if (safeIfSumSmaller(costOf(newSubset), costOf(newComplement), best.cost)) {
                best.set(dataOf(newSubset), dataOf(newComplement));
            }
        }
    }
 
    void computePartialSolution(NodeArray<DWMSplit>& split, node v, SubsetEnumerator<node>& subset,
            const List<node>& terminals) {
        List<node> newSubset;
        makeKey(newSubset, v);
 
        if (!m_map.member(newSubset)) { // not already defined
            T oldCost = costOf(terminals);
            DWMData best;
            auto addPair = [&](node v1, node v2, T dist) {
                best.add(NodePair(v1, v2), dist);
                if (m_pred[v1][v2] == nullptr) {
                    best.invalidate();
                }
            };
            for (node w : m_G.nodes) {
                T dist = m_distance[v][w];
                if (m_terminalSubset.hasMember(w)) {
                    // we attach edge vw to tree containing terminal w
                    if (safeIfSumSmaller(oldCost, dist, best.cost)) {
                        best.clear();
                        best.add(dataOf(terminals));
                        addPair(v, w, dist);
                    }
                } else {
                    // we attach edge vw to tree split[w]
                    OGDF_ASSERT(!m_terminalSubset.hasMember(v));
                    computeSplit(split, w, subset);
                    if (safeIfSumSmaller(split[w].cost, dist, best.cost)) {
                        best.clear();
                        best.add(split[w].subgraph1);
                        best.add(split[w].subgraph2);
                        if (v != w) {
                            addPair(v, w, dist);
                        }
                    }
                }
            }
            m_map.fastInsert(newSubset, best);
        }
    }
 
    template<typename CONTAINER>
    void computePartialSolutions(const CONTAINER& nodeContainer) {
        List<node> terminals;
        m_terminalSubset.list(terminals);
        SubsetEnumerator<node> subset(terminals); // done here because of linear running time
        NodeArray<DWMSplit> split(m_G);
        for (node v : nodeContainer) {
            if (!m_terminalSubset.hasMember(v)) {
                computePartialSolution(split, v, subset, terminals);
            }
        }
    }
 
    void initializeMap() {
        for (node v : m_G.nodes) {
            for (node t : m_terminals) {
                if (t != v) {
                    List<node> key;
                    key.pushBack(t);
                    if (v->index() < t->index()) {
                        key.pushFront(v);
                    } else {
                        key.pushBack(v);
                    }
 
                    if (!m_map.member(key)) { // not already defined
                        if (m_pred[t][v] == nullptr) {
                            m_map.fastInsert(key, DWMData(m_distance[t][v]));
                            OGDF_ASSERT(!dataOf(key)->valid() || m_distance[t][v] == 0);
                        } else {
                            NodePairs nodepairs;
                            nodepairs.push(NodePair(key.front(), key.back()));
                            m_map.fastInsert(key, DWMData(m_distance[t][v], nodepairs));
                            OGDF_ASSERT(dataOf(key)->valid());
                        }
                    }
                }
            }
        }
    }
 
    T getSteinerTreeFor(const DWMData& data, EdgeWeightedGraphCopy<T>& tree) const {
        T cost(0);
        if (data.valid()) {
            // add edges
            for (auto nodepair : data.nodepairs) {
                node uO = nodepair.source;
                node vO = nodepair.target;
                node uC = tree.copy(uO);
                node vC = tree.copy(vO);
                if (uC == nullptr) {
                    uC = tree.newNode(uO);
                }
                if (vC == nullptr) {
                    vC = tree.newNode(vO);
                }
#ifdef OGDF_FULL_COMPONENT_GENERATION_TERMINAL_SSSP_AWARE
                const T dist = m_isTerminal[uO] ? m_distance[uO][vO] : m_distance[vO][uO];
#else
                const T dist = m_distance[uO][vO];
#endif
                tree.newEdge(uC, vC, dist);
                cost += dist;
            }
 
            // recurse
            for (const DWMData* subgraph : data.subgraphs) {
                cost += getSteinerTreeFor(*subgraph, tree);
            }
        }
 
        OGDF_ASSERT(isAcyclicUndirected(tree));
        return cost;
    }
 
public:
    FullComponentGeneratorDreyfusWagner(const EdgeWeightedGraph<T>& G, const List<node>& terminals,
            const NodeArray<bool>& isTerminal, const NodeArray<NodeArray<T>>& distance,
            const NodeArray<NodeArray<edge>>& pred)
        : m_G(G)
        , m_terminals(terminals)
        , m_isTerminal(isTerminal)
        , m_distance(distance)
        , m_pred(pred)
        , m_terminalSubset(m_terminals)
        , m_map(1 << 22) // we initially allocate 4MB*sizeof(DWMData) for hashing
    {
        initializeMap();
    }
 
    void call(int restricted) {
        OGDF_ASSERT(restricted >= 2);
        Math::updateMin(restricted, m_terminals.size());
        for (m_terminalSubset.begin(2, restricted - 1); m_terminalSubset.valid();
                m_terminalSubset.next()) {
            if (m_terminalSubset.size() != restricted - 1) {
                computePartialSolutions(m_G.nodes);
            } else { // maximal terminal subset
                // save time by only adding terminals instead of all nodes
                computePartialSolutions(m_terminals);
            }
        }
    }
 
    T getSteinerTreeFor(const List<node>& terminals, EdgeWeightedGraphCopy<T>& tree) const {
        tree.setOriginalGraph(m_G);
        T cost(getSteinerTreeFor(*dataOf(terminals), tree));
        OGDF_ASSERT(isTree(tree));
        return cost;
    }
 
    bool isValidComponent(const EdgeWeightedGraphCopy<T>& graph) const {
        if (graph.empty()) {
            return false;
        }
        for (node v : graph.nodes) {
            if (m_isTerminal[graph.original(v)] // is a terminal
                    && v->degree() > 1) { // but not a leaf
                return false;
            }
        }
        return true;
    }
};
 
template<typename T>
class FullComponentGeneratorDreyfusWagner<T>::SortedNodeListHashFunc {
    static const unsigned int c_prime = 0x7fffffff; // mersenne prime 2**31 - 1
    // would be nicer: 0x1fffffffffffffff; // mersenne prime 2**61 - 1
    const int m_random;
 
public:
    SortedNodeListHashFunc() : m_random(randomNumber(2, c_prime - 1)) { }
 
    unsigned int hash(const List<node>& key) const {
        unsigned int hash = 0;
        for (node v : key) {
            hash = (hash * m_random + v->index()) % c_prime;
        }
        return hash;
    }
};
 
}
}