1/** \file
2 * \brief Implementation of an approxmiation algorithm for
3 * Steiner tree problems provided by Michel X. Goemans and
4 * David P. Williamson.
5 *
6 * \author Tilo Wiedera
7 *
8 * \par License:
9 * This file is part of the Open Graph Drawing Framework (OGDF).
10 *
11 * \par
12 * Copyright (C)<br>
13 * See README.md in the OGDF root directory for details.
14 *
15 * \par
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * Version 2 or 3 as published by the Free Software Foundation;
19 * see the file LICENSE.txt included in the packaging of this file
20 * for details.
21 *
22 * \par
23 * This program is distributed in the hope that it will be useful,
24 * but WITHOUT ANY WARRANTY; without even the implied warranty of
25 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
26 * GNU General Public License for more details.
27 *
28 * \par
29 * You should have received a copy of the GNU General Public
30 * License along with this program; if not, see
31 * http://www.gnu.org/copyleft/gpl.html
32 */
33
34#pragma once
35
36#include <limits>
37
38#include <ogdf/graphalg/steiner_tree/EdgeWeightedGraphCopy.h>
39#include <ogdf/module/MinSteinerTreeModule.h>
40#include <ogdf/basic/DisjointSets.h>
41
42//#define OGDF_MINSTEINERTREE_PRIMAL_DUAL_LOGGING
43
44namespace ogdf {
45
46/**
47 * @brief Primal-Dual approximation algorithm for Steiner tree problems.
48 * Yields a guaranteed approximation factor of two.
49 *
50 * This algorithm was first described by Michel X. Goemans and David P. Williamson in
51 * "A General Approximation Technique for Constrained Forest Problems",
52 * SIAM Journal on Computing, 24:296-317, 1995.
53 *
54 * @ingroup ga-steiner
55 */
56template<typename T>
57class MinSteinerTreePrimalDual : public MinSteinerTreeModule<T> {
58private:
59 const EdgeWeightedGraph<T> *m_pGraph;
60 const List<node> *m_pTerminals;
61 const NodeArray<bool> *m_pIsTerminal;
62 const T MAX_VALUE = std::numeric_limits<T>::max();
63
64 NodeArray<int> m_componentMapping;
65 DisjointSets<> *m_pComponents;
66 HashArray<int, ListIterator<int>> m_activeComponentIterators;
67 List<int> m_activeComponents;
68 double m_lowerBound;
69 NodeArray<double> m_priorities;
70
71 /**
72 * Merges two disjoint components
73 *
74 * @param v representative of the first component
75 * @param w representative of the second component
76 */
77 void mergeComponents(const node v, const node w);
78
79 /**
80 * Marks the specified component as active.
81 *
82 * @param component the component to be activated.
83 * @return
84 */
85 void makeActive(int component);
86
87 /**
88 * Returns whether the given component is active.
89 *
90 * @return true if the component is active, false otherwise
91 */
92 bool isActive(int component) const;
93
94 /**
95 * Finds the biggest set including node v.
96 *
97 * @param v the representative of the set to find
98 */
99 int getComponent(const node v) const;
100
101 /**
102 * Idendifies the next edge with a tight-to-be packing constraint.
103 *
104 * @param nextEdge the found edge
105 * @return the adjusted weight (aka epsilon) for the found edge
106 */
107 double getNextEdge(edge *nextEdge);
108
109 /**
110 * Must be called after merging any two components.
111 * Will update all the priorities of all active edges by epsilon.
112 *
113 * @param eps the value of the last tight edge
114 */
115 void updatePriorities(double eps);
116
117 /**
118 * Initializes all required datastructes.
119 */
120 void init();
121
122protected:
123 /**
124 * Builds a minimum Steiner tree given a weighted graph and a list of terminals
125 *
126 * \param G
127 * The weighted input graph
128 * \param terminals
129 * The list of terminal nodes
130 * \param isTerminal
131 * A bool array of terminals
132 * \param finalSteinerTree
133 * The final Steiner tree
134 *
135 * \return
136 * The objective value (sum of edge costs) of the final Steiner tree
137 */
138 virtual T computeSteinerTree(const EdgeWeightedGraph<T> &G, const List<node> &terminals, const NodeArray<bool> &isTerminal,
139 EdgeWeightedGraphCopy<T> *&finalSteinerTree) override;
140
141public:
142 virtual T call(const EdgeWeightedGraph<T> &G, const List<node> &terminals, const NodeArray<bool> &isTerminal, EdgeWeightedGraphCopy<T> *&finalSteinerTree) override
143 {
144 m_lowerBound = 0;
145 return MinSteinerTreeModule<T>::call(G, terminals, isTerminal, finalSteinerTree);
146 }
147
148 /**
149 * Returns the lower bound calculated while
150 * solving the last problem. Will return 0 if no problem
151 * was solved before.
152 */
153 double getLastLowerBound() const;
154};
155
156template<typename T>
157void MinSteinerTreePrimalDual<T>::init()
158{
159 m_activeComponentIterators.clear();
160 m_activeComponents.clear();
161 m_componentMapping.init(*m_pGraph);
162 m_priorities.init(*m_pGraph, 0);
163}
164
165template<typename T>
166int MinSteinerTreePrimalDual<T>::getComponent(const node v) const
167{
168 return m_pComponents->find(m_componentMapping[v]);
169}
170
171template<typename T>
172bool MinSteinerTreePrimalDual<T>::isActive(int component) const
173{
174 return m_activeComponentIterators[component].valid();
175}
176
177template<typename T>
178void MinSteinerTreePrimalDual<T>::makeActive(int comp)
179{
180 m_activeComponentIterators[comp] = m_activeComponents.pushBack(comp);
181}
182
183template<typename T>
184void MinSteinerTreePrimalDual<T>::mergeComponents(const node v, const node w)
185{
186
187 int compV = getComponent(v);
188 int compW = getComponent(w);
189
190 // remove former components
191 if(isActive(compV)) {
192 m_activeComponents.del(m_activeComponentIterators[compV]);
193 }
194 if(isActive(compW)) {
195 m_activeComponents.del(m_activeComponentIterators[compW]);
196 }
197
198 // craete new component
199 int compNew = m_pComponents->link(compV, compW);
200 if(!m_activeComponents.empty()) {
201 makeActive(compNew);
202 }
203}
204
205template<typename T>
206void MinSteinerTreePrimalDual<T>::updatePriorities(double eps)
207{
208 List<node> nodes;
209 m_pGraph->allNodes(nodes);
210 for(node v : nodes) {
211 if(isActive(getComponent(v))) {
212 m_priorities(v) += eps;
213 }
214 }
215}
216
217template<typename T>
218double MinSteinerTreePrimalDual<T>::getNextEdge(edge *nextEdge)
219{
220 double result = MAX_VALUE;
221 *nextEdge = nullptr;
222
223 List<edge> edges;
224 m_pGraph->allEdges(edges);
225 for(edge e : edges) {
226 node v = e->source();
227 node w = e->target();
228 int compV = getComponent(v);
229 int compW = getComponent(w);
230 if(compV != compW) { // spanning different components ?
231 double value = m_pGraph->weight(e) - m_priorities(v) - m_priorities(w);
232 int divisor = isActive(compV) + isActive(compW);
233 if(divisor == 0) {
234 value = MAX_VALUE;
235 } else {
236 value /= divisor;
237 }
238 if(*nextEdge == nullptr || value < result) {
239 *nextEdge = e;
240 result = value;
241 }
242 }
243 }
244 return result;
245}
246
247template<typename T>
248T MinSteinerTreePrimalDual<T>::computeSteinerTree(const EdgeWeightedGraph<T> &G, const List<node> &terminals, const NodeArray<bool> &isTerminal, EdgeWeightedGraphCopy<T> *&finalSteinerTree)
249{
250 m_pGraph = &G;
251 m_pTerminals = &terminals;
252 m_pIsTerminal = &isTerminal;
253 DisjointSets<> components;
254 m_pComponents = &components;
255
256 finalSteinerTree = new EdgeWeightedGraphCopy<T>();
257 finalSteinerTree->createEmpty(*m_pGraph);
258
259 init();
260
261 // initialize components
262 List<node> nodes;
263 m_pGraph->allNodes(nodes);
264 for(node v : nodes) {
265 int comp = m_pComponents->makeSet();
266 m_componentMapping[v] = comp;
267 if(isTerminal(v)) {
268 makeActive(comp);
269 }
270 }
271
272#ifdef OGDF_MINSTEINERTREE_PRIMAL_DUAL_LOGGING
273 std::cout << "Goemans primal-dual starting..." << std::endl;
274 std::cout << "terminals:";
275 for(node v : *m_pTerminals) {
276 std::cout << " " << v;
277 }
278 std::cout << std::endl;
279
280 std::cout << "loop starting... " << std::endl;
281#endif
282
283 T result = 0;
284 while(!m_activeComponents.empty()) {
285#ifdef OGDF_MINSTEINERTREE_PRIMAL_DUAL_LOGGING
286 std::cout << "active component exists" << std::endl;
287#endif
288 // idendify next edge
289 edge minEdge = nullptr;
290 double minValue = getNextEdge(&minEdge);
291 OGDF_ASSERT(minEdge != nullptr);
292
293#ifdef OGDF_MINSTEINERTREE_PRIMAL_DUAL_LOGGING
294 std::cout << "minEdge found: " << minEdge << ", weight is " << m_pGraph->weight(minEdge) << ", adjusted weight is " << minValue << std::endl;
295#endif
296 node v = minEdge->source();
297 node w = minEdge->target();
298
299 // include nodes in Steiner tree
300 if(finalSteinerTree->copy(v) == nullptr) {
301 finalSteinerTree->newNode(v);
302 }
303 if(finalSteinerTree->copy(w) == nullptr) {
304 finalSteinerTree->newNode(w);
305 }
306
307 // include edge in Steiner tree
308 T weight = m_pGraph->weight(minEdge);
309 result += weight;
310 finalSteinerTree->newEdge(minEdge, weight);
311
312 m_lowerBound += m_activeComponents.size() * minValue;
313
314 mergeComponents(v, w);
315
316 updatePriorities(minValue);
317 }
318 result -= this->pruneAllDanglingSteinerPaths(*finalSteinerTree, *m_pIsTerminal);
319
320#ifdef OGDF_MINSTEINERTREE_PRIMAL_DUAL_LOGGING
321 std::cout << "calculation finished!" << std::endl;
322#endif
323 return result;
324}
325
326template<typename T>
327double MinSteinerTreePrimalDual<T>::getLastLowerBound() const
328{
329 return m_lowerBound;
330}
331
332}
333