Open
Graph Drawing
Framework

 v. 2023.09 (Elderberry)
 

Loading...
Searching...
No Matches
Approximation.h
Go to the documentation of this file.
1
32#pragma once
33
38
39#include <memory>
40
41//#define OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
42
43namespace ogdf {
44namespace steiner_tree {
45namespace goemans {
46
48template<typename T>
54
55 const double m_eps;
56
57 std::minstd_rand m_rng;
58
60 struct TemporaryEdges : ArrayBuffer<edge> {
63
65 edge add(node v, node w, T cost, int capacity) {
66 edge e = m_blowupGraph.newEdge(v, w, cost, capacity);
67 this->push(e);
68 return e;
69 }
70
72 ~TemporaryEdges() { m_blowupGraph.delEdges(*this); }
73
74 private:
76 };
77
81 return maxFlow.computeValue(blowupGraph.capacities(), blowupGraph.getSource(),
82 blowupGraph.getTarget());
83 }
84
87 int findComponentAndMaxBasis(std::unique_ptr<ArrayBuffer<std::pair<node, int>>>& maxBasis,
89 // there should always be saturated flow to the component roots
90 // (contracted matroid)
91 EdgeArray<int> lB(blowupGraph.getGraph(), 0);
92 for (adjEntry adj : blowupGraph.getSource()->adjEntries) {
93 const edge e = adj->theEdge();
94 lB[e] = blowupGraph.getCapacity(e);
95 }
96
97 // compute weights of core edges and add source->core edges
99 EdgeArray<double> cost(blowupGraph.getGraph(), 0);
100 for (node v : blowupGraph.core()) {
101 // compute weight of core edge
102 double weight = blowupGraph.computeCoreWeight(v);
103
104 // add edges from source to core edges v
105 edge e = sourceCoreEdges.add(blowupGraph.getSource(), v, 0,
106 blowupGraph.getCoreCapacity(v));
107 cost[e] = -weight;
108 }
109
110 NodeArray<int> supply(blowupGraph.getGraph(), 0);
111 EdgeArray<int> flow(blowupGraph.getGraph());
113
114 for (int id = 1; id <= gamma.size(); ++id) {
115 /* We want to find maximum-weight basis B \in B^K_Q
116 * B_Q = minimal edge set to remove after contracting Q to obtain feasible solution (blowup graph)
117 * = bases of gammoid M_Q
118 * B^K_Q = { B \in B_Q | B \subseteq K } [K is the set of core edges]
119 * = bases of gammoid M^K_Q
120 * M^K_Q is gammoid "obtained by restricting M_Q to K"
121 * M_Q = M'_Q / X' (M'_Q contracted by X') is a gammoid
122 * M'_Q = gammoid from X \cup X' to Y in D' with arc-disjointness instead of vertex-disjointness
123 * D' is D like for separation but without source node and with interim node v_e for each edge e
124 * X = inserted nodes v_e in each edge in blowup graph
125 * X' = the (arbitrary) roots of each component
126 * Y = Q \cup {t}
127 * that means:
128 * I (subset of X) is an independent set of M_Q
129 * <=> (if X' is always an independent set in M'_Q ... which we assume here)
130 * I \cup X' is an independent set of M'_Q
131 * Restricting to K means:
132 * only put v_e in X with e \in K (that is, we only need to *generate* v_e if e \in K)
133 * Hence:
134 * - we generate D'^K (this is blowupGraph)
135 * - compute the max flow from X^K \cup X' to Q \cup {t}
136 * - ASSERT that X' is "saturated"!
137 * - check which subset of X is saturated -> these are the nodes representing the edge set we need
138 */
139 // add edges from component's terminals to target
141 for (node t : gamma.terminals(id)) {
142 Q_to_target.add(t, blowupGraph.getTarget(), 0,
143 blowupGraph.getLCM() * blowupGraph.getY());
144 // the last value is an upper bound for the capacity
145 }
146 // TODO: we could also use static edges from all terminals to the target
147 // and just change their capacities each time; needs extra data structures
148
149 std::unique_ptr<ArrayBuffer<std::pair<node, int>>> basis(
150 new ArrayBuffer<std::pair<node, int>>);
151
152 int rank = gammoidGetRank(blowupGraph);
153 OGDF_ASSERT(rank >= blowupGraph.getY() + blowupGraph.getLCM());
154 supply[blowupGraph.getSource()] = rank;
155 supply[blowupGraph.getTarget()] = -rank;
156
157#ifdef OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
158 std::cout << "Computing min-cost flow for blowup component " << id << " of "
159 << gamma.size() << std::endl;
160 std::cout << " * terminals of component are " << gamma.terminals(id) << std::endl;
161 for (node v : blowupGraph.getGraph().nodes) {
162 if (supply[v] > 0) {
163 std::cout << " * supply node " << v << " with supply " << supply[v] << std::endl;
164 }
165 if (supply[v] < 0) {
166 std::cout << " * demand node " << v << " with demand " << -supply[v] << std::endl;
167 }
168 }
169 for (edge e : blowupGraph.getGraph().edges) {
170 std::cout << " * edge " << e << " with cost " << blowupGraph.getCost(e)
171 << " and flow bounds [" << lB[e] << ", " << blowupGraph.getCapacity(e)
172 << "]" << std::endl;
173 }
174#endif
175
176 // find maximum weight basis
177#ifdef OGDF_DEBUG
178 bool feasible =
179#endif
180 mcf.call(blowupGraph.getGraph(), lB, blowupGraph.capacities(), cost, supply,
181 flow);
182 OGDF_ASSERT(feasible);
183 OGDF_ASSERT(mcf.checkComputedFlow(blowupGraph.getGraph(), lB, blowupGraph.capacities(),
184 cost, supply, flow));
185
186 double weight(0);
187 for (edge e : sourceCoreEdges) {
188 if (flow[e] > 0) {
189 basis->push(std::make_pair(e->target(), flow[e]));
190 weight -= flow[e] * cost[e];
191 }
192 }
193
194#ifdef OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
195 std::cout << "Basis weight is " << weight << std::endl
196 << "Checking if " << gamma.cost(id) << "(component cost) * "
197 << blowupGraph.getLCM() << "(lcm) <= " << weight << "(basis weight)"
198 << std::endl;
199#endif
200
201 // we choose the component with max cost*N <= weight
202 if (gamma.cost(id) * blowupGraph.getLCM() <= weight + m_eps) {
203#ifdef OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
204 std::cout << "Using basis because it is feasible." << std::endl;
205#endif
206 maxBasis.swap(basis);
207 return id;
208 }
209 }
210 return 0; // no component chosen, fail
211 }
212
216 std::unique_ptr<ArrayBuffer<std::pair<node, int>>>& maxBasis,
217 const BlowupGraph<T>& blowupGraph, const BlowupComponents<T>& gamma) {
218 // find cheapest component
219 int compId = 0;
220 double cost = 0;
221 for (int id = 1; id <= gamma.size(); ++id) {
222 if (gamma.cost(id) > cost) {
223 cost = gamma.cost(id);
224 compId = id;
225 }
226 }
227 // use all core edges as basis
228 maxBasis.reset(new ArrayBuffer<std::pair<node, int>>());
229 for (node v : blowupGraph.core()) {
230 maxBasis->push(std::make_pair(v, blowupGraph.getCoreCapacity(v)));
231 }
232 return compId;
233 }
234
237 const edge rootEdge) {
238 OGDF_ASSERT(blowupGraph.isTerminal(rootEdge->source()));
240 stack.push(rootEdge->target());
241 while (!stack.empty()) {
242 const node v = stack.popRet();
243 if (blowupGraph.isTerminal(v)) {
244 continue;
245 }
246 const node vO = blowupGraph.getOriginal(v);
247 if (vO) {
248 isNewTerminal[vO] = true;
249 }
250 for (adjEntry adj : v->adjEntries) {
251 const node w = adj->theEdge()->target();
252 if (v != w) { // outgoing edge
253 stack.push(w);
254 }
255 }
256 }
257 }
258
262#ifdef OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
263 std::cout << "Remove basis from blowup graph" << std::endl;
264#endif
265 // remove B from K (K := K \ B) and from blowup graph (X := X - B)
266 // and, while at it, remove cleanup edges from blowup graph (X := X - F)
267 // and fix components that have no incoming edges
268 ArrayBuffer<Prioritized<node, int>> fractionalCoreEdges; // we defer fractional basis elements
269 for (auto p : basis) {
270 const node v = p.first;
271 const int count = p.second;
272 int origCap = blowupGraph.getCoreCapacity(v);
274#ifdef OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
275 std::cout << " * node " << v << " with count " << count << " (of " << origCap << ")"
276 << std::endl;
277#endif
278 if (count < origCap) { // only remove a fraction?
280#ifdef OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
281 std::cout << " -> deferred because fractional" << std::endl;
282#endif
283 } else {
284 // we are deleting the core edge from the whole component
285 blowupGraph.delCore(v);
286 blowupGraph.removeBasis(v);
287#ifdef OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
288 std::cout << " -> done" << std::endl;
289#endif
290 }
291 }
292 fractionalCoreEdges.quicksort(); // sort decreasing by flow value
293#ifdef OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
294 if (!fractionalCoreEdges.empty()) {
295 std::cout << "Deferred core edges:" << std::endl;
296 }
297#endif
298 for (auto p : fractionalCoreEdges) {
299 const node v = p.item();
300 const int count = -p.priority();
301 int origCap = blowupGraph.getCoreCapacity(v);
302#ifdef OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
303 std::cout << " * node " << v << " with count " << count << " of " << origCap << std::endl;
304#endif
306 // copy (split) the component
307 blowupGraph.copyComponent(blowupGraph.findRootEdge(v), count, origCap - count);
308 // we are deleting the core edge from the whole component
309 blowupGraph.delCore(v);
310 blowupGraph.removeBasis(v);
311 }
312 }
313
314public:
317 const NodeArray<bool>& isTerminal,
319 const std::minstd_rand& rng, double eps = 1e-8)
320 : m_G(G)
321 , m_isTerminal(isTerminal)
322 , m_terminals(terminals)
324 , m_eps(eps)
325 , m_rng(rng) { }
326
330};
331
332template<typename T>
334#ifdef OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
335 std::cout << "Start solving based on LP solution" << std::endl;
336 int iteration = 0;
337#endif
339 BlowupGraph<T> blowupGraph(m_G, m_terminals, m_fullCompStore, cer, m_eps);
340
341 while (blowupGraph.terminals().size() > 1) { // T is not a Steiner tree
342#ifdef OGDF_STEINER_TREE_GOEMANS_APPROXIMATION_LOGGING
343 std::cout << "Iteration " << ++iteration << " with " << blowupGraph.terminals().size()
344 << " terminals" << std::endl;
345#endif
346 BlowupComponents<T> gamma(blowupGraph); // Gamma(X)
347
349
350 // take a component Q in Gamma(X)
351 std::unique_ptr<ArrayBuffer<std::pair<node, int>>> maxBasis;
352 int compId = blowupGraph.getY() > 0
353 ? findComponentAndMaxBasis(maxBasis, blowupGraph, gamma)
354 : findCheapestComponentAndRemainingBasis(maxBasis, blowupGraph, gamma);
355 OGDF_ASSERT(compId != 0);
356
357 // add component Q to T
358 addComponent(isNewTerminal, blowupGraph, gamma.rootEdge(compId));
359
360 // remove (maybe fractional) basis and do all the small things necessary for update
361 removeFractionalBasisAndCleanup(*maxBasis, blowupGraph, gamma);
362
363 // contract (X := X / Q)
364 blowupGraph.contract(gamma.terminals(compId));
365
366 if (blowupGraph.terminals().size() > 1) {
367 blowupGraph.updateSpecialCapacities();
368 }
369 }
370}
371
372}
373}
374}
Definition of the ogdf::steiner_tree::goemans::BlowupComponents class template.
Definition of ogdf::steiner_tree::goemans::CoreEdgeRandomSpanningTree class template.
Declaration and implementation of Goldberg-Tarjan max-flow algorithm with global relabeling and gap r...
Definition of ogdf::MinCostFlowReinelt class template.
Class for adjacency list elements.
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(edge e)
Puts a new element in the buffer.
int capacity() const
Returns the current capacity of the datastructure. Note that this value is rather irrelevant if the a...
Dynamic arrays indexed with edges.
Definition EdgeArray.h:125
Class for the representation of edges.
Definition Graph_d.h:300
node target() const
Returns the target node of the edge.
Definition Graph_d.h:338
node source() const
Returns the source node of the edge.
Definition Graph_d.h:335
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).
static bool checkComputedFlow(const Graph &G, const EdgeArray< int > &lowerBound, const EdgeArray< int > &upperBound, const EdgeArray< TCost > &cost, const NodeArray< int > &supply, const EdgeArray< int > &flow, TCost &value)
checks if a computed flow is a feasible solution to the given problem instance.
Computes a min-cost flow using a network simplex method.
virtual bool call(const Graph &G, const EdgeArray< int > &lowerBound, const EdgeArray< int > &upperBound, const EdgeArray< TCost > &cost, const NodeArray< int > &supply, EdgeArray< int > &flow, NodeArray< TCost > &dual) override
Computes a min-cost flow in the directed graph G using a network simplex method.
Dynamic arrays indexed with nodes.
Definition NodeArray.h:125
Class for the representation of nodes.
Definition Graph_d.h:177
internal::GraphObjectContainer< AdjElement > adjEntries
The container containing all entries in the adjacency list of this node.
Definition Graph_d.h:208
Augments any data elements of type X with keys of type Priority. This class is also its own Comparer.
Definition comparer.h:291
A data structure to store full components with extra data for each component.
The actual 1.39-approximation algorithm by Goemans et al. with a set of terminalized nodes as result.
int findCheapestComponentAndRemainingBasis(std::unique_ptr< ArrayBuffer< std::pair< node, int > > > &maxBasis, const BlowupGraph< T > &blowupGraph, const BlowupComponents< T > &gamma)
For the end of the algorithm: find cheapest component and choose all remaining core edges as basis.
void removeFractionalBasisAndCleanup(ArrayBuffer< std::pair< node, int > > &basis, BlowupGraph< T > &blowupGraph, BlowupComponents< T > &gamma)
Remove a given basis and cleanup, the basis may be given fractionally.
const EdgeWeightedGraph< T > & m_G
const FullComponentWithExtraStore< T, double > & m_fullCompStore
all enumerated full components, with solution
int gammoidGetRank(const BlowupGraph< T > &blowupGraph) const
Computes the rank of the gammoid (given by the blowup graph)
Approximation(const EdgeWeightedGraph< T > &G, const List< node > &terminals, const NodeArray< bool > &isTerminal, const FullComponentWithExtraStore< T, double > &fullCompStore, const std::minstd_rand &rng, double eps=1e-8)
Initialize everything.
void addComponent(NodeArray< bool > &isNewTerminal, const BlowupGraph< T > &blowupGraph, const edge rootEdge)
Add a component of the blowup graph to the final solution (by changing nonterminals to terminals)
const double m_eps
epsilon for double operations
int findComponentAndMaxBasis(std::unique_ptr< ArrayBuffer< std::pair< node, int > > > &maxBasis, BlowupGraph< T > &blowupGraph, const BlowupComponents< T > &gamma)
Finds the best component and its maximum-weight basis in the given blowup graph with given core and w...
void solve(NodeArray< bool > &isNewTerminal)
Perform the actual approximation algorithm on the LP solution.
Obtain and provides information about components in a given blowup graph.
A special-purpose blowup graph for gammoid computation: directed, with special source and target,...
Definition BlowupGraph.h:50
bool isLoopFree(const Graph &G)
Returns true iff G contains no self-loop.
#define OGDF_ASSERT(expr)
Assert condition expr. See doc/build.md for more information.
Definition basic.h:41
static MultilevelBuilder * getDoubleFactoredZeroAdjustedMerger()
The namespace for all OGDF objects.
Add edges into a blowup graph and delete them on destruction.
edge add(node v, node w, T cost, int capacity)
Add a temporary edge to the blowup graph.
TemporaryEdges(BlowupGraph< T > &blowupGraph)
Construct object for a specific blowupGraph.