1/** \file
2 * \brief Bandit test suite for Steiner tree algorithms
3 *
4 * \author Tilo Wiedera
5 *
6 * \par License:
7 * This file is part of the Open Graph Drawing Framework (OGDF).
8 *
9 * \par
10 * Copyright (C)<br>
11 * See README.md in the OGDF root directory for details.
12 *
13 * \par
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * Version 2 or 3 as published by the Free Software Foundation;
17 * see the file LICENSE.txt included in the packaging of this file
18 * for details.
19 *
20 * \par
21 * This program is distributed in the hope that it will be useful,
22 * but WITHOUT ANY WARRANTY; without even the implied warranty of
23 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 * GNU General Public License for more details.
25 *
26 * \par
27 * You should have received a copy of the GNU General Public
28 * License along with this program; if not, see
29 * http://www.gnu.org/copyleft/gpl.html
30 */
31
32#include <string>
33#include <vector>
34#include <ogdf/fileformats/GraphIO.h>
35#include <ogdf/graphalg/MinSteinerTreeDirectedCut.h>
36#include <ogdf/graphalg/MaxFlowEdmondsKarp.h>
37#include <ogdf/graphalg/MinSteinerTreeKou.h>
38#include <ogdf/graphalg/MinSteinerTreeMehlhorn.h>
39#include <ogdf/graphalg/MinSteinerTreeRZLoss.h>
40#include <ogdf/graphalg/MinSteinerTreeZelikovsky.h>
41#include <ogdf/graphalg/MinSteinerTreeShore.h>
42#include <ogdf/graphalg/MinSteinerTreePrimalDual.h>
43#include <ogdf/graphalg/MinSteinerTreeDualAscent.h>
44#include <ogdf/graphalg/MinSteinerTreeGoemans139.h>
45#include <resources.h>
46
47template<typename T>
48struct ModuleData {
49 //! a human-readable name/description of the module
50 std::string name;
51 //! the Steiner tree module to be tested
52 std::unique_ptr<MinSteinerTreeModule<T>> alg;
53 //! the approximation factor of this algorithm, needed for validating the results
54 double ratio;
55 //! the sizes (number of nodes) of the random graphs to test
56 std::vector<int> sizes;
57};
58
59template<typename T>
60using Modules = std::vector<ModuleData<T>>;
61
62template<typename T>
63static void addModule(Modules<T>& modules, const std::string& name, MinSteinerTreeModule<T>* alg, double ratio, std::vector<int> sizes = {35, 50}) {
64 modules.emplace_back(ModuleData<T>{name, std::unique_ptr<MinSteinerTreeModule<T>>(alg), ratio, sizes});
65}
66
67/**
68 * Generates a new graph with an optimal Steiner tree.
69 * Only very basic graphs are generated
70 * to guarantee the optimality of the resulting Steiner tree.
71 *
72 * \param n
73 * number of nodes
74 * \param graph
75 * the resulting graph
76 * \param terminals
77 * this list will hold all terminals
78 * \param isTerminal
79 * stores which node is a terminal
80 * \param tree
81 * an optimal Steiner tree for this graph.
82 */
83template<typename T>
84T randomOptimalSteiner(
85 int n,
86 EdgeWeightedGraph<T> &graph,
87 List<node> &terminals,
88 NodeArray<bool> &isTerminal,
89 EdgeWeightedGraphCopy<T> &tree
90)
91{
92 OGDF_ASSERT(n >= 4);
93
94 terminals.clear();
95
96 int numberOfTerminals = max(3, randomNumber(n/4, n/2));
97 int numberOfNonterminals = n - numberOfTerminals;
98 int numberOfEdges = randomNumber(numberOfTerminals-1 + numberOfNonterminals*2, (n*(n-1))/2);
99
100 randomTree(graph, numberOfTerminals);
101 isTerminal.init(graph, false);
102 for (node v : graph.nodes) {
103 if (v->degree() == 1) {
104 isTerminal[v] = true;
105 }
106 }
107 for (edge e : graph.edges) {
108 graph.setWeight(e, 1);
109 }
110
111 tree.init(graph);
112 T result = tree.numberOfEdges();
113
114 for(int i = numberOfTerminals-1; i < numberOfEdges;) {
115 node v = graph.chooseNode();
116 node u = graph.chooseNode([&](node w) { return w != v; });
117 OGDF_ASSERT(u != nullptr);
118
119 if(numberOfNonterminals > 0) {
120 node w = graph.newNode();
121 graph.newEdge(v, w, n);
122 graph.newEdge(w, u, n);
123 numberOfNonterminals--;
124 i += 2;
125 }
126 else {
127 if (graph.searchEdge(v, u) == nullptr
128 && graph.searchEdge(u, v) == nullptr) {
129 graph.newEdge(v, u, n);
130 i++;
131 }
132 }
133 }
134
135 MinSteinerTreeModule<T>::getTerminals(terminals, graph, isTerminal);
136
137 OGDF_ASSERT(terminals.size() <= numberOfTerminals);
138 OGDF_ASSERT(graph.numberOfEdges() == numberOfEdges);
139 OGDF_ASSERT(tree.numberOfNodes() == numberOfTerminals);
140 OGDF_ASSERT(tree.numberOfEdges() == numberOfTerminals - 1);
141 OGDF_ASSERT(tree.weight(tree.firstEdge()) == 1);
142 OGDF_ASSERT(tree.weight(tree.lastEdge()) == 1);
143 OGDF_ASSERT(graph.numberOfNodes() == n);
144 OGDF_ASSERT(isSimpleUndirected(graph));
145 OGDF_ASSERT(isConnected(graph));
146 return result;
147}
148
149/**
150 * Generates a random Steiner tree instance.
151 *
152 * \param n number of nodes
153 * \param graph the resulting graph
154 * \param terminals this list will hold all terminals
155 * \param isTerminal stores which node is a terminal
156 */
157template<typename T>
158void randomSteinerTreeInstance(
159 int n,
160 EdgeWeightedGraph<T> &graph,
161 List<node> &terminals,
162 NodeArray<bool> &isTerminal)
163{
164 OGDF_ASSERT(n >= 3);
165
166 randomSimpleConnectedGraph(graph, n, randomNumber(2*n - 3, n*(n-1)/2));
167 int numberOfTerminals = max(3, randomNumber(n/4, 2*n/3));
168
169 for (edge e : graph.edges) {
170 graph.setWeight(e, randomNumber(1, 100));
171 }
172
173 Array<node> nodes;
174 graph.allNodes(nodes);
175 nodes.permute();
176
177 terminals.clear();
178 isTerminal.init(graph, false);
179 for (int i = 0; i < numberOfTerminals; ++i) {
180 const node v = nodes[i];
181 isTerminal[v] = true;
182 terminals.pushBack(v);
183 }
184}
185
186/**
187 * Test if modules generates a valid/reasonable Steiner tree for a graph with given number of nodes
188 */
189template<typename T>
190void testModuleOnRandomGraph(MinSteinerTreeModule<T> &alg, int n, double factor = 0)
191{
192 it("generates a valid Steiner tree for a random graph of " + to_string(n) + " nodes", [&]() {
193 EdgeWeightedGraph<T> graph;
194 NodeArray<bool> isTerminal(graph, false);
195 List<node> terminals;
196
197 randomSteinerTreeInstance(n, graph, terminals, isTerminal);
198 std::cout << " (" << terminals.size() << " terminals, " << graph.numberOfEdges() << " edges)";
199
200 EdgeWeightedGraphCopy<T>* make_solution;
201 T returnedCost = alg.call(graph, terminals, isTerminal, make_solution);
202
203 std::unique_ptr<EdgeWeightedGraphCopy<T>> solution(make_solution);
204
205 for (node v : solution->nodes) {
206 AssertThat(solution->original(v), !IsNull());
207 }
208
209 T actualCost(0);
210 for (edge e : solution->edges) {
211 AssertThat(solution->original(e), !IsNull());
212 actualCost += solution->weight(e);
213 }
214
215 AssertThat(actualCost, Equals(returnedCost));
216 AssertThat(MinSteinerTreeModule<T>::isSteinerTree(graph, terminals, isTerminal, *solution), Equals(true));
217 });
218
219 it("finds a reasonable Steiner tree for a graph of " + to_string(n) + " nodes", [&]() {
220 EdgeWeightedGraph<T> graph;
221 EdgeWeightedGraphCopy<T> tree;
222 NodeArray<bool> isTerminal(graph, false);
223 List<node> terminals;
224
225 T cost = randomOptimalSteiner<T>(n, graph, terminals, isTerminal, tree);
226 std::cout << " (" << terminals.size() << " terminals, " << graph.numberOfEdges() << " edges)";
227
228 EdgeWeightedGraphCopy<T>* make_solution;
229 T algCost = alg.call(graph, terminals, isTerminal, make_solution);
230 std::unique_ptr<EdgeWeightedGraphCopy<T>> solution(make_solution);
231
232 AssertThat(MinSteinerTreeModule<T>::isSteinerTree(graph, terminals, isTerminal, *solution), IsTrue());
233
234 // only check optimum approximation
235 // for algorithms with factor of 2 or better
236 if(factor >= 1 && factor <= 2) {
237 AssertThat(algCost, Equals(cost));
238 AssertThat(solution->numberOfNodes(), Equals(tree.numberOfNodes()));
239 AssertThat(solution->numberOfEdges(), Equals(tree.numberOfEdges()));
240
241 List<node> nodes;
242 tree.allNodes(nodes);
243 for(node v : nodes) {
244 AssertThat(solution->copy(tree.original(v)), !IsNull());
245 }
246
247 List<edge> edges;
248 tree.allEdges(edges);
249 for(edge e : edges) {
250 AssertThat(solution->copy(tree.original(e)), !IsNull());
251 }
252 }
253 });
254}
255
256/**
257 * Tests one subclass of MinSteinerTreeModule for a specific type.
258 */
259template<typename T>
260void testModule(const ModuleData<T>& module)
261{
262 describe(module.name, [&]() {
263 for (int n : module.sizes) {
264 testModuleOnRandomGraph(*module.alg, n, module.ratio);
265 }
266
267 for_each_file("steiner", [&](const ResourceFile* file){
268 // optimal solution value is extracted from the filename
269 string filename = file->name();
270 string tmp = filename.substr(0, filename.length() - 4);
271 T opt(0);
272 auto pos = tmp.find_last_of('.');
273 if (pos != tmp.npos) {
274 tmp = tmp.substr(tmp.find_last_of('.') + 1);
275 std::stringstream ss(tmp);
276 ss >> opt;
277 }
278
279 it("yields correct results on " + file->fullPath() + " (optimum is " + (opt == 0 ? "unknown" : to_string(opt)) + ")", [&] {
280 EdgeWeightedGraph<T> graph;
281 List<node> terminals;
282 NodeArray<bool> isTerminal;
283
284 std::stringstream is{file->data()};
285 GraphIO::readSTP(graph, terminals, isTerminal, is);
286
287 EdgeWeightedGraphCopy<T>* make_solution;
288 T algCost = module.alg->call(graph, terminals, isTerminal, make_solution);
289 std::unique_ptr<EdgeWeightedGraphCopy<T>> solution(make_solution);
290
291 AssertThat(MinSteinerTreeModule<T>::isSteinerTree(graph, terminals, isTerminal, *solution), IsTrue());
292 if (opt > 0) {
293 AssertThat(algCost, IsGreaterThan(opt) || Equals(opt));
294 if (module.ratio != 0) {
295 AssertThat(algCost, IsLessThan(module.ratio*opt) || Equals(module.ratio*opt));
296 }
297 }
298 });
299 });
300 });
301}
302
303struct MaxFlowFactoryBase {
304 virtual MaxFlowModule<double>* create() = 0;
305 virtual ~MaxFlowFactoryBase() = default;
306};
307template<typename MaxFlowModuleType>
308struct MaxFlowFactory : MaxFlowFactoryBase {
309 MaxFlowModule<double>* create() override {
310 return new MaxFlowModuleType();
311 }
312};
313
314/**
315 * Registers one instance of the MinSteinerTreeDirectedCut class for each of its variants
316 */
317template <typename T>
318static void
319registerDirectedCutVariants(Modules<T> &modules)
320{
321 using AlgPair = std::pair<MaxFlowFactoryBase*, std::string>;
322 std::unique_ptr<MaxFlowFactoryBase> flowEK(new MaxFlowFactory<MaxFlowEdmondsKarp<double>>());
323 std::unique_ptr<MaxFlowFactoryBase> flowGT(new MaxFlowFactory<MaxFlowEdmondsKarp<double>>());
324
325 using BoolPair = std::pair<bool, std::string>;
326 struct VerboseTrueFalse : public std::vector<BoolPair> {
327 VerboseTrueFalse(std::string&& trueString, std::string&& falseString = "") : std::vector<BoolPair>({{true, trueString}, {false, falseString}}) {}
328 };
329
330 for (auto maxFlow : {AlgPair{flowEK.get(), "Edmonds-Karp"}, AlgPair{flowGT.get(), "Goldberg-Tarjan"}}) {
331 for (auto useBackCuts : VerboseTrueFalse{", back cuts"}) {
332 for (auto useMinCardinalityCuts : VerboseTrueFalse{", min cardinality cuts"}) {
333 for (auto useNestedCuts : VerboseTrueFalse{", nested cuts"}) {
334 for (auto useExtraConstraints : VerboseTrueFalse{"all extra constraints", "only necessary constraints"}) {
335 MinSteinerTreeDirectedCut<T> *alg = new MinSteinerTreeDirectedCut<T>();
336
337 std::stringstream ss;
338 ss << "DirectedCut";
339
340 alg->setMaxFlowModule(maxFlow.first->create());
341 ss << ", " << maxFlow.second;
342
343 alg->useBackCuts(useBackCuts.first);
344 ss << useBackCuts.second;
345
346 alg->useMinCardinalityCuts(useMinCardinalityCuts.first);
347 ss << useMinCardinalityCuts.second;
348
349 alg->useNestedCuts(useNestedCuts.first);
350 ss << useNestedCuts.second;
351
352 alg->useDegreeConstraints(useExtraConstraints.first);
353 alg->useFlowBalanceConstraints(useExtraConstraints.first);
354 alg->useGSEC2Constraints(useExtraConstraints.first);
355 alg->useIndegreeEdgeConstraints(useExtraConstraints.first);
356 ss << ", " << useExtraConstraints.second;
357
358 addModule(modules, ss.str(), alg, 1, {12, 30});
359 }
360 }
361 }
362 }
363 }
364}
365
366/**
367 * Registers one instance of the MinSteinerTreeZelikovsky class for each of its variants
368 */
369template <typename T>
370static void
371registerZelikovskyVariants(Modules<T> &modules)
372{
373 using WCalc = std::tuple<std::string, typename MinSteinerTreeZelikovsky<T>::WinCalculation>;
374 using TGen = std::tuple<std::string, typename MinSteinerTreeZelikovsky<T>::TripleGeneration>;
375 using TRed = std::tuple<std::string, typename MinSteinerTreeZelikovsky<T>::TripleReduction>;
376 using SCalc = std::tuple<std::string, typename MinSteinerTreeZelikovsky<T>::SaveCalculation>;
377 using Pass = std::tuple<std::string, typename MinSteinerTreeZelikovsky<T>::Pass>;
378 using APSPStrategy = std::tuple<std::string, bool>;
379
380 std::vector<WCalc> winCalculations = {
381 WCalc("absolute win function", MinSteinerTreeZelikovsky<T>::WinCalculation::absolute),
382 WCalc("relative win function", MinSteinerTreeZelikovsky<T>::WinCalculation::relative)
383 };
384 std::vector<TGen> tripleGenStrategies = {
385 TGen("exhaustive triple generation", MinSteinerTreeZelikovsky<T>::TripleGeneration::exhaustive),
386 TGen("Voronoi triple generation", MinSteinerTreeZelikovsky<T>::TripleGeneration::voronoi),
387 TGen("direct triple generation", MinSteinerTreeZelikovsky<T>::TripleGeneration::ondemand)
388 };
389 std::vector<TRed> tripleReductStrategies = {
390 TRed("enabled reduction", MinSteinerTreeZelikovsky<T>::TripleReduction::on),
391 TRed("disabled reduction", MinSteinerTreeZelikovsky<T>::TripleReduction::off),
392 };
393 std::vector<SCalc> saveCalculations = {
394 SCalc("static enumeration save calculation", MinSteinerTreeZelikovsky<T>::SaveCalculation::staticEnum),
395 SCalc("static LCATree save calculation", MinSteinerTreeZelikovsky<T>::SaveCalculation::staticLCATree),
396 SCalc("dynamic LCATree save calculation", MinSteinerTreeZelikovsky<T>::SaveCalculation::dynamicLCATree),
397 SCalc("hybrid save calculation", MinSteinerTreeZelikovsky<T>::SaveCalculation::hybrid)
398 };
399 std::vector<Pass> passes = {
400 Pass("one-pass", MinSteinerTreeZelikovsky<T>::Pass::one),
401 Pass("multi-pass", MinSteinerTreeZelikovsky<T>::Pass::multi)
402 };
403 std::vector<APSPStrategy> apspStrategies = {
404 APSPStrategy("forced APSP", true),
405 APSPStrategy("SSSP", false),
406 };
407
408 std::vector<typename decltype(winCalculations)::size_type>
409 choice = { 0, 0, 0, 0, 0, 0 },
410 maxchoice = {
411 winCalculations.size(),
412 tripleGenStrategies.size(),
413 tripleReductStrategies.size(),
414 saveCalculations.size(),
415 passes.size(),
416 apspStrategies.size(),
417 };
418 enum indices {
419 winIdx = 0,
420 tgenIdx,
421 tredIdx,
422 saveIdx,
423 passIdx,
424 apspIdx,
425 };
426
427 auto nextChoice = [&]() {
428 bool overflow;
429 unsigned int i = 0;
430 do {
431 overflow = false;
432 ++choice[i];
433 choice[i] %= maxchoice[i];
434 if (choice[i] == 0) {
435 ++i;
436 overflow = true;
437 }
438 } while (overflow && i < choice.size());
439 return !overflow;
440 };
441
442 do {
443 string desc = "Zelikovsky: ";
444
445 MinSteinerTreeZelikovsky<T> *module = new MinSteinerTreeZelikovsky<T>();
446
447 Pass pass = passes[choice[passIdx]];
448 desc += std::get<0>(pass);
449 module->pass(std::get<1>(pass));
450
451 SCalc saveCalc = saveCalculations[choice[saveIdx]];
452 desc += ", " + std::get<0>(saveCalc);
453 module->saveCalculation(std::get<1>(saveCalc));
454
455 TGen tripleGen = tripleGenStrategies[choice[tgenIdx]];
456 desc += ", " + std::get<0>(tripleGen);
457 module->tripleGeneration(std::get<1>(tripleGen));
458
459 TRed tripleRed = tripleReductStrategies[choice[tredIdx]];
460 desc += ", " + std::get<0>(tripleRed);
461 module->tripleReduction(std::get<1>(tripleRed));
462
463 WCalc winCalc = winCalculations[choice[winIdx]];
464 desc += ", " + std::get<0>(winCalc);
465 module->winCalculation(std::get<1>(winCalc));
466
467 APSPStrategy apspStrategy = apspStrategies[choice[apspIdx]];
468 desc += ", " + std::get<0>(apspStrategy);
469 module->forceAPSP(std::get<1>(apspStrategy));
470
471 // check for invalid configurations
472 if (module->tripleGeneration() == MinSteinerTreeZelikovsky<T>::TripleGeneration::ondemand
473 && ((module->winCalculation() != MinSteinerTreeZelikovsky<T>::WinCalculation::absolute)
474 || (module->saveCalculation() == MinSteinerTreeZelikovsky<T>::SaveCalculation::hybrid)
475 || (module->tripleReduction() == MinSteinerTreeZelikovsky<T>::TripleReduction::off)
476 || (module->pass() == MinSteinerTreeZelikovsky<T>::Pass::one))) {
477 delete module;
478 } else {
479 addModule(modules, desc, module, 11/6.0);
480 };
481 } while (nextChoice());
482}
483
484/**
485 * Registers one instance of the MinSteinerTreeRZLoss class for each of its variants
486 */
487template <typename T>
488static void
489registerRZLossVariants(Modules<T> &modules)
490{
491 // RZLoss for different maximum component sizes
492 for(int i = 2; i < 6; i++) {
493 MinSteinerTreeRZLoss<T> *alg = new MinSteinerTreeRZLoss<T>();
494 int maxCompSize = i;
495 std::string info = "";
496 // APSP is only being used for maximum component size of 3
497 if(i == 2) {
498 alg->forceAPSP(true);
499 info = " and forced APSP";
500 maxCompSize = 3;
501 }
502 alg->setMaxComponentSize(maxCompSize);
503 addModule(modules, "RZLoss with maximum component size of " + to_string(maxCompSize) + info, alg, 2, {14, 25});
504 }
505}
506
507/**
508 * Registers one instance of the MinSteinerTreeGoemans139 class for each of its variants
509 */
510template <typename T>
511static void
512registerGoemans139Variants(Modules<T> &modules)
513{
514 // Goemans139 for different maximum component sizes
515 for(int i = 2; i < 6; i++) {
516 // and for standard and stronger LP relaxation
517 for (int strongerLP = 0; strongerLP < 2; ++strongerLP) {
518 for (int use2approx = 0; use2approx < 2; ++use2approx) {
519 MinSteinerTreeGoemans139<T> *alg = new MinSteinerTreeGoemans139<T>();
520 int maxCompSize = i;
521 std::string info = "Goemans139 with maximum component size ";
522 if(i == 2) {
523 alg->forceAPSP();
524 maxCompSize = 3;
525 info += "3 (enforced APSP)";
526 } else {
527 info += to_string(maxCompSize);
528 }
529 alg->setMaxComponentSize(maxCompSize);
530 if (strongerLP) {
531 alg->separateCycles();
532 info += " using stronger LP";
533 }
534 if (use2approx) {
535 alg->use2Approximation();
536 info += " with upper bound";
537 }
538 addModule(modules, info, alg, 2, {14, 25});
539 }
540 }
541 }
542}
543
544/**
545 * Registers a complete Steiner test suite for a given
546 * template parameter, like int or double.
547 */
548template<typename T>
549void registerSuite(const std::string typeName)
550{
551 describe("for graphs with " + typeName + "-typed costs:", [] {
552 Modules<T> modules;
553 addModule(modules, "DirectedCut default", new MinSteinerTreeDirectedCut<T>(), 1);
554 addModule(modules, "Kou", new MinSteinerTreeKou<T>(), 2);
555 addModule(modules, "Mehlhorn", new MinSteinerTreeMehlhorn<T>(), 2);
556 addModule(modules, "RZLoss default", new MinSteinerTreeRZLoss<T>(), 2);
557 addModule(modules, "Goemans139 default", new MinSteinerTreeGoemans139<T>(), 2);
558 addModule(modules, "Takahashi", new MinSteinerTreeTakahashi<T>(), 2);
559 addModule(modules, "Shore", new MinSteinerTreeShore<T>(), 1, {10, 20});
560 addModule(modules, "Primal-Dual", new MinSteinerTreePrimalDual<T>(), 2);
561 addModule(modules, "DualAscent", new MinSteinerTreeDualAscent<T>(), 0);
562 addModule(modules, "Zelikovsky default", new MinSteinerTreeZelikovsky<T>(), 11/6.0);
563
564 registerDirectedCutVariants<T>(modules);
565 registerZelikovskyVariants<T>(modules);
566 registerRZLossVariants<T>(modules);
567 registerGoemans139Variants<T>(modules);
568
569 // register suites
570 for (auto& module : modules) {
571 testModule<T>(module);
572 module.alg.reset();
573 }
574 });
575}
576
577go_bandit([](){
578 describe("Steiner tree algorithms", []() {
579 registerSuite<int>("int");
580 registerSuite<double>("double");
581 });
582});
583