MinSteinerTreeZelikovsky.h

Go to the documentation of this file.

 
#pragma once
 
#include <ogdf/basic/ArrayBuffer.h>
#include <ogdf/basic/Graph.h>
#include <ogdf/basic/List.h>
#include <ogdf/basic/basic.h>
#include <ogdf/basic/comparer.h>
#include <ogdf/basic/extended_graph_alg.h>
#include <ogdf/graphalg/MinSteinerTreeModule.h>
#include <ogdf/graphalg/steiner_tree/EdgeWeightedGraphCopy.h>
#include <ogdf/graphalg/steiner_tree/Full3ComponentGeneratorEnumeration.h>
#include <ogdf/graphalg/steiner_tree/Full3ComponentGeneratorVoronoi.h>
#include <ogdf/graphalg/steiner_tree/Save.h>
#include <ogdf/graphalg/steiner_tree/SaveDynamic.h>
#include <ogdf/graphalg/steiner_tree/SaveEnum.h>
#include <ogdf/graphalg/steiner_tree/SaveStatic.h>
#include <ogdf/graphalg/steiner_tree/Triple.h>
#include <ogdf/graphalg/steiner_tree/common_algorithms.h>
 
#include <limits>
 
namespace ogdf::steiner_tree {
template<typename T>
class Full3ComponentGeneratorModule;
} // namespace ogdf::steiner_tree
 
namespace ogdf {
template<typename T>
class EdgeWeightedGraph;
 
template<typename T>
class MinSteinerTreeZelikovsky : public MinSteinerTreeModule<T> {
public:
    template<typename TYPE>
    using Save = steiner_tree::Save<TYPE>;
    template<typename TYPE>
    using Triple = steiner_tree::Triple<TYPE>;
 
    enum class WinCalculation {
        absolute, 
        relative 
    };
 
    enum class TripleGeneration {
        exhaustive, 
        voronoi, 
        ondemand 
    };
 
    enum class TripleReduction {
        on, 
        off 
    };
 
    enum class SaveCalculation {
        staticEnum,
        staticLCATree,
        dynamicLCATree,
        hybrid
    };
 
    enum class Pass {
        one,
        multi
    };
 
    MinSteinerTreeZelikovsky(WinCalculation wc = WinCalculation::absolute,
            TripleGeneration tg = TripleGeneration::voronoi,
            SaveCalculation sc = SaveCalculation::hybrid, TripleReduction tr = TripleReduction::on,
            Pass pass = Pass::multi)
        : m_winCalculation(wc)
        , m_tripleGeneration(tg)
        , m_saveCalculation(sc)
        , m_tripleReduction(tr)
        , m_pass(pass)
        , m_ssspDistances(true) { }
 
    virtual ~MinSteinerTreeZelikovsky() { }
 
    virtual T call(const EdgeWeightedGraph<T>& G, const List<node>& terminals,
            const NodeArray<bool>& isTerminal, EdgeWeightedGraphCopy<T>*& finalSteinerTree) override {
        m_triplesGenerated = 0;
        m_tripleLookUps = 0;
        m_triplesContracted = 0;
        return MinSteinerTreeModule<T>::call(G, terminals, isTerminal, finalSteinerTree);
    }
 
    void forceAPSP(bool force = true) { m_ssspDistances = !force; }
 
    void winCalculation(WinCalculation wc) { m_winCalculation = wc; }
 
    WinCalculation winCalculation() const { return m_winCalculation; }
 
    void tripleGeneration(TripleGeneration tg) { m_tripleGeneration = tg; }
 
    TripleGeneration tripleGeneration() const { return m_tripleGeneration; }
 
    void tripleReduction(TripleReduction tr) { m_tripleReduction = tr; }
 
    TripleReduction tripleReduction() const { return m_tripleReduction; }
 
    void saveCalculation(SaveCalculation sv) { m_saveCalculation = sv; }
 
    SaveCalculation saveCalculation() const { return m_saveCalculation; }
 
    void pass(Pass p) { m_pass = p; }
 
    Pass pass() const { return m_pass; }
 
    long numberOfGeneratedTriples() const { return m_triplesGenerated; }
 
    long numberOfContractedTriples() const { return m_triplesContracted; }
 
    long numberOfTripleLookUps() const { return m_tripleLookUps; }
 
protected:
    virtual T computeSteinerTree(const EdgeWeightedGraph<T>& G, const List<node>& terminals,
            const NodeArray<bool>& isTerminal, EdgeWeightedGraphCopy<T>*& finalSteinerTree) override;
 
    void computeDistanceMatrix();
 
    inline void generateTriple(node u, node v, node w, node center, const T& minCost,
            const Save<T>& save) {
        const double gain = save.gain(u, v, w);
        const double win = calcWin(gain, minCost);
        if (tripleReduction() == TripleReduction::off || win > 0) {
            ++m_triplesGenerated;
            OGDF_ASSERT(center);
            Triple<T> triple(u, v, w, center, minCost, win);
            m_triples.pushBack(triple);
        }
    }
 
    inline void generateTriples(const Save<T>& save,
            const steiner_tree::Full3ComponentGeneratorModule<T>& fcg) {
        fcg.call(*m_originalGraph, *m_terminals, *m_isTerminal, m_distance, m_pred,
                [this, &save](node u, node v, node w, node center, T minCost) {
                    generateTriple(u, v, w, center, minCost, save);
                });
    }
 
    inline void generateTriples(const Save<T>& save) {
        OGDF_ASSERT(tripleGeneration() != TripleGeneration::ondemand);
        if (tripleGeneration() == TripleGeneration::voronoi) {
            steiner_tree::Full3ComponentGeneratorVoronoi<T> fcg;
            generateTriples(save, fcg);
        } else {
            OGDF_ASSERT(tripleGeneration() == TripleGeneration::exhaustive);
            steiner_tree::Full3ComponentGeneratorEnumeration<T> fcg;
            generateTriples(save, fcg);
        }
    }
 
    void contractTriple(const Triple<T>& triple, Save<T>& save, NodeArray<bool>& isNewTerminal) {
        ++m_triplesContracted;
        save.update(triple);
        isNewTerminal[triple.z()] = true;
    }
 
    void tripleOnDemand(Save<T>& save, NodeArray<bool>& isNewTerminal);
 
    bool findBestTripleForCenter(node center, const Save<T>& save, Triple<T>& maxTriple) const {
        bool updated = false; // return value
 
        // find s0, nearest terminal to center
        T best = std::numeric_limits<T>::max();
        node s0 = nullptr;
        for (node s : *m_terminals) {
            T tmp = m_distance[s][center];
            if (best > tmp) {
                best = tmp;
                s0 = s;
            }
        }
        OGDF_ASSERT(s0);
        OGDF_ASSERT(m_pred[s0][center]);
 
        // find s1 maximizing save(s0, s1) - d(center, s1)
        node s1 = nullptr;
        T save1Val(0);
        for (node s : *m_terminals) {
            if (s != s0 && m_pred[s][center] != nullptr) {
                OGDF_ASSERT(m_distance[s][center] != std::numeric_limits<T>::max());
                T tmpVal = save.saveWeight(s, s0);
                T tmp = tmpVal - m_distance[s][center];
                if (!s1 || best < tmp) {
                    best = tmp;
                    s1 = s;
                    save1Val = tmpVal;
                }
            }
        }
        if (s1) {
            OGDF_ASSERT(m_pred[s1][center]);
            node s2 = nullptr;
            T save2Val(0);
            const edge save1 = save.saveEdge(s0, s1);
            for (node s : *m_terminals) {
                if (s != s0 && s != s1 && m_pred[s][center] != nullptr) {
                    OGDF_ASSERT(m_distance[s][center] != std::numeric_limits<T>::max());
                    const edge tmp = save.saveEdge(s0, s);
                    save2Val = save.saveWeight(tmp == save1 ? s1 : s0, s);
                    T tmpWin = save1Val + save2Val - m_distance[s0][center] - m_distance[s1][center]
                            - m_distance[s][center];
                    if (!s2 || best < tmpWin) {
                        best = tmpWin;
                        s2 = s;
                    }
                }
            }
 
            if (s2 // it may happen that s2 does not exist
                    && best > maxTriple.win()) { // best win is better than previous best; also positive
                OGDF_ASSERT(m_pred[s2][center]);
                maxTriple.s0(s0);
                maxTriple.s1(s1);
                maxTriple.s2(s2);
                maxTriple.z(center);
                maxTriple.win(best);
                //maxTriple.cost(save1Val + save2Val - win); not needed
                updated = true;
            }
        }
        return updated;
    }
 
    void multiPass(Save<T>& save, NodeArray<bool>& isNewTerminal);
 
    void onePass(Save<T>& save, NodeArray<bool>& isNewTerminal);
 
    double calcWin(double gain, T cost) const {
        switch (winCalculation()) {
        case WinCalculation::relative:
            return gain / cost - 1.0;
        case WinCalculation::absolute:
        default:
            return gain - cost;
        }
    }
 
    void generateInitialTerminalSpanningTree(EdgeWeightedGraphCopy<T>& steinerTree) {
        // generate complete graph
        for (node v : *m_terminals) {
            steinerTree.newNode(v);
        }
        for (node u : steinerTree.nodes) {
            const NodeArray<T>& dist = m_distance[steinerTree.original(u)];
            for (node v = u->succ(); v; v = v->succ()) {
                steinerTree.newEdge(u, v, dist[steinerTree.original(v)]);
            }
        }
        // compute MST
        makeMinimumSpanningTree(steinerTree, steinerTree.edgeWeights());
    }
 
private:
    WinCalculation m_winCalculation; 
    TripleGeneration m_tripleGeneration; 
    SaveCalculation m_saveCalculation; 
    TripleReduction m_tripleReduction; 
    Pass m_pass; 
    bool m_ssspDistances; 
 
    const EdgeWeightedGraph<T>* m_originalGraph; 
    const NodeArray<bool>* m_isTerminal; 
    const List<node>* m_terminals; 
    NodeArray<NodeArray<T>> m_distance; 
    NodeArray<NodeArray<edge>> m_pred; 
    List<Triple<T>> m_triples; 
 
    long m_triplesGenerated; 
    long m_triplesContracted; 
    long m_tripleLookUps; 
};
 
template<typename T>
T MinSteinerTreeZelikovsky<T>::computeSteinerTree(const EdgeWeightedGraph<T>& G,
        const List<node>& terminals, const NodeArray<bool>& isTerminal,
        EdgeWeightedGraphCopy<T>*& finalSteinerTree) {
    OGDF_ASSERT(tripleGeneration()
                    != TripleGeneration::ondemand // combinations that only work with ondemand:
            || (winCalculation() == WinCalculation::absolute
                    && saveCalculation() != SaveCalculation::hybrid && pass() != Pass::one));
 
    m_originalGraph = &G;
    m_terminals = &terminals;
    m_isTerminal = &isTerminal;
 
    // build the complete terminal graph and a distance matrix
    NodeArray<bool> isNewTerminal(G, false);
    for (node v : *m_terminals) {
        isNewTerminal[v] = true;
    }
 
    if (terminals.size() >= 3) {
        computeDistanceMatrix();
 
        // init terminal-spanning tree and its save-edge data structure
        EdgeWeightedGraphCopy<T> steinerTree; // the terminal-spanning tree to be modified
        steinerTree.setOriginalGraph(G);
        generateInitialTerminalSpanningTree(steinerTree);
 
        Save<T>* save = nullptr;
        switch (saveCalculation()) {
        case SaveCalculation::staticEnum:
            save = new steiner_tree::SaveEnum<T>(steinerTree);
            break;
        case SaveCalculation::staticLCATree:
            save = new steiner_tree::SaveStatic<T>(steinerTree);
            break;
        case SaveCalculation::dynamicLCATree:
        case SaveCalculation::hybrid:
            save = new steiner_tree::SaveDynamic<T>(steinerTree);
            break;
        }
        OGDF_ASSERT(save);
 
        if (tripleGeneration() == TripleGeneration::ondemand) { // ondemand triple generation
            tripleOnDemand(*save, isNewTerminal);
        } else { // exhaustive or voronoi
            // triple generation phase
            if (saveCalculation() == SaveCalculation::hybrid) {
                steiner_tree::SaveEnum<T> saveTriple(steinerTree);
                generateTriples(saveTriple);
            } else {
                generateTriples(*save);
            }
            // contraction phase
            if (pass() == Pass::multi) {
                multiPass(*save, isNewTerminal);
            } else { // pass() == one
                onePass(*save, isNewTerminal);
            }
        }
        delete save;
 
        // cleanup
        m_triples.clear();
    }
 
    // obtain final Steiner Tree using (MST-based) Steiner tree approximation algorithm
    return steiner_tree::obtainFinalSteinerTree(G, isNewTerminal, isTerminal, finalSteinerTree);
}
 
template<typename T>
void MinSteinerTreeZelikovsky<T>::computeDistanceMatrix() {
    if (m_ssspDistances) {
        MinSteinerTreeModule<T>::allTerminalShortestPaths(*m_originalGraph, *m_terminals,
                *m_isTerminal, m_distance, m_pred);
    } else {
        MinSteinerTreeModule<T>::allPairShortestPaths(*m_originalGraph, *m_isTerminal, m_distance,
                m_pred);
    }
}
 
template<typename T>
void MinSteinerTreeZelikovsky<T>::tripleOnDemand(Save<T>& save, NodeArray<bool>& isNewTerminal) {
    Triple<T> maxTriple;
    ArrayBuffer<node> nonterminals;
    MinSteinerTreeModule<T>::getNonterminals(nonterminals, *m_originalGraph, *m_isTerminal);
    do {
        maxTriple.win(0);
        for (node center : nonterminals) {
            if (findBestTripleForCenter(center, save, maxTriple)) {
                ++m_triplesGenerated;
            }
        }
 
        if (maxTriple.win() > 0) {
            contractTriple(maxTriple, save, isNewTerminal);
        }
    } while (maxTriple.win() > 0);
}
 
template<typename T>
void MinSteinerTreeZelikovsky<T>::onePass(Save<T>& save, NodeArray<bool>& isNewTerminal) {
    m_triples.quicksort(GenericComparer<Triple<T>, double>(
            [](const Triple<T>& x) -> double { return -x.win(); }));
 
    for (const Triple<T>& t : m_triples) {
        ++m_tripleLookUps;
        if (calcWin(double(save.gain(t.s0(), t.s1(), t.s2())), t.cost()) > 0) {
            contractTriple(t, save, isNewTerminal);
        }
    }
}
 
template<typename T>
void MinSteinerTreeZelikovsky<T>::multiPass(Save<T>& save, NodeArray<bool>& isNewTerminal) {
    double win = 0;
    ListIterator<Triple<T>> maxTriple;
 
    do {
        win = 0;
        ListIterator<Triple<T>> nextIt;
        for (ListIterator<Triple<T>> it = m_triples.begin(); it.valid(); it = nextIt) {
            nextIt = it.succ();
            ++m_tripleLookUps;
            Triple<T>& t = *it;
            t.win(calcWin(double(save.gain(t.s0(), t.s1(), t.s2())), t.cost()));
            if (t.win() > win) {
                win = t.win();
                maxTriple = it;
            } else {
                if (tripleReduction() == TripleReduction::on && t.win() <= 0) {
                    m_triples.del(it);
                }
            }
        }
 
        if (win > 0) {
            contractTriple(*maxTriple, save, isNewTerminal);
            m_triples.del(maxTriple);
        }
    } while (win > 0);
}
 
}