# OpenGraph DrawingFramework

v. 2023.09 (Elderberry)

Searching...
No Matches
LPRelaxationSER.h
Go to the documentation of this file.
1
32#pragma once
33
34#include <coin/CoinPackedMatrix.hpp>
35
37#include <ogdf/external/coin.h>
39
40
41//#define OGDF_STEINERTREE_LPRELAXATIONSER_LOGGING
42//#define OGDF_STEINERTREE_LPRELAXATIONSER_OUTPUT_LP
43#define OGDF_STEINERTREE_LPRELAXATIONSER_SEPARATE_CONNECTED_COMPONENTS // this is faster
44#define OGDF_STEINERTREE_LPRELAXATIONSER_SEPARATE_YVAR_CONSTRAINTS // if not defined: generate all yvar constraints in the beginning
45
46namespace ogdf {
47namespace steiner_tree {
48
51template<typename T>
57
60 double* m_objective;
63
64 const T m_upperBound;
66#ifdef OGDF_STEINERTREE_LPRELAXATIONSER_SEPARATE_CONNECTED_COMPONENTS
68#endif
69
70 const double m_eps;
71
73 void generateProblem();
80
84
88
92
95 bool separate();
96
100 node& target, GraphCopy& G, EdgeArray<double>& capacity, int& cutsFound) {
101 G.createEmpty(m_G);
102 capacity.init(G);
103 source = G.newNode();
104 for (node t : m_terminals) { // generate all terminals
105 G.newNode(t);
106 }
107 target = G.newNode();
108 for (int j = 0; j < activeComponents.size(); ++j) {
109 const int i = activeComponents[j];
110 const double cap = m_fullCompStore.extra(i);
111 const Array<node>& terminals = m_fullCompStore.terminals(i);
112 // take the first terminal as root
113 // XXX: we may generate parallel edges but it's ok
114 const auto it0 = terminals.begin();
115 const node rC = G.copy(*it0);
116 capacity[G.newEdge(source, rC)] = cap;
117 if (terminals.size() > 2) {
118 const node v = G.newNode();
119 capacity[G.newEdge(rC, v)] = cap;
120 for (auto it = it0 + 1; it != terminals.end(); ++it) {
121 const node w = G.copy(*it);
122 capacity[G.newEdge(v, w)] = cap;
123 }
124 } else { // exactly two terminals: we do not need the inner Steiner node
125 capacity[G.newEdge(rC, G.copy(*terminals.rbegin()))] = cap;
126 }
127 }
128 double y_R = 0;
129 // TODO: perhaps better to compute y_v before
130 // add edges to target and compute y_R
131 for (node t : m_terminals) {
132 const node v = G.copy(t);
133 OGDF_ASSERT(v);
134 // compute y_v, simply the sum of all x_C where C contains v and then - 1
135 double y_v(-1);
137 if (adj->twinNode() != source) {
139 }
140 }
141
142#ifdef OGDF_STEINERTREE_LPRELAXATIONSER_SEPARATE_YVAR_CONSTRAINTS
143 if (y_v < -m_eps) {
145 ++cutsFound;
146 } else
147#endif
148 if (y_v > 0) {
149 capacity[G.newEdge(v, target)] = y_v;
150 y_R += y_v;
151 }
152 }
153#if 0
154 // just for output of blow-up graph
155 for (edge e : G.edges) {
156 if (G.original(e->source())) {
157 std::cout << " T:" << G.original(e->source());
158 } else {
159 std::cout << " " << e->source();
160 }
161 std::cout << " -> ";
162 if (G.original(e->target())) {
163 std::cout << "T:" << G.original(e->target());
164 } else {
165 std::cout << e->target();
166 }
167 std::cout << " " << capacity[e] << "\n";
168 }
169#endif
170
171 return y_R;
172 }
173
174public:
187 T upperBound = 0, int cliqueSize = 0, double eps = 1e-8)
188 : m_G(G)
189 , m_isTerminal(isTerminal)
190 , m_terminals(terminals)
192 , m_osiSolver(CoinManager::createCorrectOsiSolverInterface())
197 , m_upperBound(upperBound)
201#endif
202 , m_eps(eps) {
205
206#ifndef OGDF_STEINERTREE_LPRELAXATIONSER_SEPARATE_YVAR_CONSTRAINTS
207 for (node t : m_terminals) {
209 }
210#endif
211 }
212
214 delete[] m_objective;
215 delete m_matrix;
216 delete[] m_lowerBounds;
217 delete[] m_upperBounds;
218 delete m_osiSolver;
219 }
220
223 bool solve();
224};
225
226template<typename T>
228 const int n = m_fullCompStore.size();
229
230 m_matrix->setDimensions(0, n);
231 for (int i = 0; i < n; ++i) {
232 m_lowerBounds[i] = 0;
233 m_upperBounds[i] = 1;
234 }
235 for (int i = 0; i < n; ++i) {
236 m_objective[i] = m_fullCompStore.cost(i);
237 }
238 m_osiSolver->loadProblem(*m_matrix, m_lowerBounds, m_upperBounds, m_objective, nullptr, nullptr);
239
240 if (m_upperBound > 0) { // add upper bound
242 row.setFull(m_fullCompStore.size(), m_objective);
244 }
245}
246
247template<typename T>
250
251 for (int i = 0; i < m_fullCompStore.size(); ++i) {
252 if (m_fullCompStore.isTerminal(i, t)) { // component spans terminal
253 row.insert(i, 1);
254 }
255 }
256
258}
259
260// add constraint if necessary and return if necessary
261template<typename T>
264 double test = 0;
265
266 for (int i = 0; i < m_fullCompStore.size(); ++i) {
267 // compute the intersection cardinality (linear time because terminal sets are sorted by index)
268 int intersectionCard = 0;
269 auto terminals = m_fullCompStore.terminals(i);
270 node* it1 = terminals.begin();
271 auto it2 = subset.begin();
272 while (it1 != terminals.end() && it2 != subset.end()) {
273 if ((*it1)->index() < (*it2)->index()) {
274 ++it1;
275 } else if ((*it1)->index() > (*it2)->index()) {
276 ++it2;
277 } else { // ==
279 ++it1;
280 ++it2;
281 }
282 }
283 // and use it as a coefficient
284 if (intersectionCard > 1) {
285 row.insert(i, intersectionCard - 1);
286 test += (intersectionCard - 1) * m_fullCompStore.extra(i);
287 }
288 }
289 if (test > subset.size() - 1) {
290 m_osiSolver->addRow(row, 0, subset.size() - 1);
291 return true;
292 }
293 return false;
294}
295
296template<typename T>
298 // we use the sum over all components C of (|C| - 1) * x_C = |R| - 1 constraint
299 // from the paper
301
302 for (int i = 0; i < m_fullCompStore.size(); ++i) {
303 row.insert(i, m_fullCompStore.terminals(i).size() - 1);
304 }
305
306 int value = m_terminals.size() - 1;
308}
309
310template<typename T>
312 bool initialIteration = true;
313
314 do {
315 if (initialIteration) {
316 m_osiSolver->initialSolve();
317#ifdef OGDF_STEINERTREE_LPRELAXATIONSER_LOGGING
318 std::cout << "Objective value " << m_osiSolver->getObjValue() << " of initial solution."
319 << std::endl;
320#endif
321 initialIteration = false;
322 } else {
323 m_osiSolver->resolve();
324#ifdef OGDF_STEINERTREE_LPRELAXATIONSER_LOGGING
325 std::cout << "Objective value " << m_osiSolver->getObjValue() << " after resolve."
326 << std::endl;
327#endif
328 }
329
330 if (!m_osiSolver->isProvenOptimal()) {
331 if (m_upperBound > 0) { // failed due to better upper bound
332 return false;
333 } else { // failed due to infeasibility
334 std::cerr << "Failed to optimize LP!" << std::endl;
335 throw(-1);
336 }
337 }
338 } while (separate());
339
340 const double* constSol = m_osiSolver->getColSolution();
341 const int numberOfColumns = m_osiSolver->getNumCols();
342
343 for (int i = 0; i < numberOfColumns; ++i) {
344 m_fullCompStore.extra(i) = constSol[i];
345 }
346
347#ifdef OGDF_STEINERTREE_LPRELAXATIONSER_OUTPUT_LP
348 m_osiSolver->writeLp(stderr);
349#endif
350
351 return true;
352}
353
354template<typename T>
356 const double* constSol = m_osiSolver->getColSolution();
357
359 for (int i = 0; i < m_fullCompStore.size(); ++i) {
360 m_fullCompStore.extra(i) = constSol[i];
361 if (m_fullCompStore.extra(i) > m_eps) {
363 }
364 }
365
366#ifdef OGDF_STEINERTREE_LPRELAXATIONSER_SEPARATE_CONNECTED_COMPONENTS
367 if (!m_separationStage) {
368 if (separateConnected(activeComponents)) {
369 return true;
370 }
371 m_separationStage = 1;
372 }
373#endif
374 if (separateMinCut(activeComponents)) {
375 return true;
376 }
377 return m_separateCliqueSize > 2 ? separateCycles(activeComponents) : false;
378}
379
380template<typename T>
382 NodeArray<int> setID(m_G, -1); // XXX: NodeArray over terminals only would be better
383 DisjointSets<> uf(m_terminals.size());
384 for (node t : m_terminals) {
385 setID[t] = uf.makeSet();
386 }
387
388 // union all nodes within one component
389 for (int j = 0; j < activeComponents.size(); ++j) {
390 auto terminals = m_fullCompStore.terminals(activeComponents[j]);
391 auto it = terminals.begin();
392 const int s1 = setID[*it];
393 for (++it; it != terminals.end(); ++it) {
395 }
396 }
397
398 if (uf.getNumberOfSets() == 1) { // solution is connected
399 return false;
400 }
401
402 Array<ArrayBuffer<node>> components(m_terminals.size());
404 for (node t : m_terminals) {
405 const int k = uf.find(setID[t]);
406 if (components[k].empty()) {
407 usedComp.push(k);
408 }
409 components[k].push(t);
410 }
411 int cutsFound = 0;
412 for (const int k : usedComp) {
414 }
415 return true;
416}
417
418template<typename T>
420 int cutsFound = 0;
421 node source;
424 EdgeArray<double> capacity;
425 double y_R = generateMinCutSeparationGraph(activeComponents, source, pseudotarget, auxG,
426 capacity, cutsFound);
427
428#ifdef OGDF_STEINERTREE_LPRELAXATIONSER_SEPARATE_YVAR_CONSTRAINTS
429 if (cutsFound > 0) {
430 return true;
431 }
432#endif
433
434 node target = auxG.newNode();
435 capacity[auxG.newEdge(pseudotarget, target)] = y_R;
436
440 for (node t : m_terminals) {
441 const node v = auxG.copy(t);
442
443 edge v_to_target = auxG.newEdge(v, target);
444 capacity[v_to_target] = std::numeric_limits<double>::max(); // XXX: smaller is better
445
446 maxFlow.init(auxG, &flow);
447
448 const double cutVal = maxFlow.computeValue(capacity, source, target);
449 if (cutVal - y_R < 1 - m_eps) {
450 minSTCut.call(auxG, capacity, flow, source, target);
452 for (node tOrig : m_terminals) {
453 const node tCopy = auxG.copy(tOrig);
454 if (tCopy && minSTCut.isInBackCut(tCopy)) {
455 subset.push(tOrig);
456 }
457 }
458
460 }
461
462 auxG.delEdge(v_to_target);
463 }
464 return cutsFound != 0;
465}
466
467template<typename T>
469 int count = 0;
470
471 // generate auxiliary graph
472 Graph G;
473 NodeArray<int> id(G);
474 for (int i : activeComponents) {
475 id[G.newNode()] = i;
476 }
477 for (node u1 : G.nodes) {
478 const int i1 = id[u1];
479 const Array<node>& terminals1 = m_fullCompStore.terminals(i1);
480 for (node u2 = u1->succ(); u2; u2 = u2->succ()) {
481 const int i2 = id[u2];
482 const Array<node>& terminals2 = m_fullCompStore.terminals(i2);
483 // compute intersection cardinality (linear time because terminal sets are sorted by index)
484 int intersectionCard = 0;
485 const node* it1 = terminals1.begin();
486 const node* it2 = terminals2.begin();
487 while (it1 != terminals1.end() && it2 != terminals2.end()) {
488 if ((*it1)->index() < (*it2)->index()) {
489 ++it1;
490 } else if ((*it1)->index() > (*it2)->index()) {
491 ++it2;
492 } else { // ==
494 if (intersectionCard == 2) {
495 G.newEdge(u1, u2);
496 break;
497 }
498 ++it1;
499 ++it2;
500 }
501 }
502 }
503 }
504
505 if (G.numberOfEdges() == 0) {
506 return false;
507 }
508
509 // now find cliques
510 Array<List<node>> degrees(m_separateCliqueSize);
511 for (node v : G.nodes) {
512 int idx = v->degree();
513 if (idx == 0) { // ignore isolated nodes
514 continue;
515 }
516 if (idx >= m_separateCliqueSize) {
517 idx = m_separateCliqueSize - 1;
518 }
519 --idx;
520 degrees[idx].pushBack(v);
521 }
522 NodeArray<bool> test(G, false);
523 for (int k = degrees.size(); k >= 2; --k) {
524 degrees[k - 2].conc(degrees[k - 1]);
525 if (degrees[k - 2].size() >= k) {
527 for (nodeSubset.begin(k); nodeSubset.valid(); nodeSubset.next()) {
528 int countEdges = (k * (k - 1)) / 2;
529 for (int j = 0; j < nodeSubset.size(); ++j) {
530 test[nodeSubset[j]] = true;
531 }
532 for (edge e : G.edges) {
533 if (test[e->source()] && test[e->target()]) {
534 countEdges -= 1;
535 }
536 }
538 if (countEdges == 0) {
539 // found clique, add constraint
540 double val(0);
542
543 for (int j = 0; j < nodeSubset.size(); ++j) {
544 int i = id[nodeSubset[j]];
545 val += m_fullCompStore.extra(i);
546 row.insert(i, 1);
547 }
548 if (val >= 1 + m_eps) {
550 ++count;
551 }
552 }
553 for (int j = 0; j < nodeSubset.size(); ++j) {
554 test[nodeSubset[j]] = false;
555 }
556 }
557 }
558 }
559 return count > 0;
560}
561
562}
563}
Definition of the FullComponentStore class template.
#define OGDF_STEINERTREE_LPRELAXATIONSER_SEPARATE_CONNECTED_COMPONENTS
Declaration of min-st-cut algorithms parameterized by max-flow alorithm.
A class that allows to enumerate k-subsets.
Definition Graph_d.h:79
An array that keeps track of the number of inserted elements; also usable as an efficient stack.
Definition ArrayBuffer.h:56
void push(E e)
Puts a new element in the buffer.
The parameterized class Array implements dynamic arrays of type E.
Definition Array.h:214
iterator begin()
Returns an iterator to the first element.
Definition Array.h:324
reverse_iterator rbegin()
Returns an reverse iterator to the last element.
Definition Array.h:342
iterator end()
Returns an iterator to one past the last element.
Definition Array.h:333
INDEX size() const
Returns the size (number of elements) of the array.
Definition Array.h:297
If you use COIN-OR, you should use this class.
Definition coin.h:43
A Union/Find data structure for maintaining disjoint sets.
Dynamic arrays indexed with edges.
Definition EdgeArray.h:125
void init()
Reinitializes the array. Associates the array with no graph.
Definition EdgeArray.h:292
Class for the representation of edges.
Definition Graph_d.h:300
Copies of graphs supporting edge splitting.
Definition GraphCopy.h:254
Data type for general directed graphs (adjacency list representation).
Definition Graph_d.h:521
Doubly linked lists (maintaining the length of the list).
Definition List.h:1435
Computes a max flow via Preflow-Push (global relabeling and gap relabeling heuristic).
Min-st-cut algorithm, that calculates the cut via maxflow.
Dynamic arrays indexed with nodes.
Definition NodeArray.h:125
Class for the representation of nodes.
Definition Graph_d.h:177
The container containing all entries in the adjacency list of this node.
Definition Graph_d.h:208
node succ() const
Returns the successor in the list of all nodes.
Definition Graph_d.h:229
Enumerator for k-subsets of a given type.
const Array< node > & terminals(int id) const
Returns the list of terminals in the full component with given id.
A data structure to store full components with extra data for each component.
Returns a reference to the extra data of this full component.
Class managing the component-based subtour elimination LP relaxation for the Steiner tree problem and...
bool separate()
Perform all available separation algorithms.
bool separateConnected(const ArrayBuffer< int > &activeComponents)
Separate to ensure that the solution is connected.
LPRelaxationSER(const EdgeWeightedGraph< T > &G, const List< node > &terminals, const NodeArray< bool > &isTerminal, FullComponentWithExtraStore< T, double > &fullCompStore, T upperBound=0, int cliqueSize=0, double eps=1e-8)
Initialize the LP.
bool solve()
Solve the LP. The solution will be written to the extra data of the full component store.
double generateMinCutSeparationGraph(const ArrayBuffer< int > &activeComponents, node &source, node &target, GraphCopy &G, EdgeArray< double > &capacity, int &cutsFound)
Generates an auxiliary multi-graph for separation (during LP solving): directed, with special source ...
bool separateCycles(const ArrayBuffer< int > &activeComponents)
Perform the separation algorithm for cycle constraints (to obtain stronger LP solution)
const double m_eps
epsilon for double operations
void generateProblem()
Generate the basic LP model.
bool addSubsetCoverConstraint(const ArrayBuffer< node > &subset)
Add subset cover constraints to the LP for a given subset of terminals.
const EdgeWeightedGraph< T > & m_G
const NodeArray< bool > & m_isTerminal
bool separateMinCut(const ArrayBuffer< int > &activeComponents)
Perform the general cut-based separation algorithm.