17 */ |
17 */ |
18 |
18 |
19 #ifndef LEMON_CYCLE_CANCELING_H |
19 #ifndef LEMON_CYCLE_CANCELING_H |
20 #define LEMON_CYCLE_CANCELING_H |
20 #define LEMON_CYCLE_CANCELING_H |
21 |
21 |
22 /// \ingroup min_cost_flow |
22 /// \ingroup min_cost_flow_algs |
23 /// |
|
24 /// \file |
23 /// \file |
25 /// \brief Cycle-canceling algorithm for finding a minimum cost flow. |
24 /// \brief Cycle-canceling algorithms for finding a minimum cost flow. |
26 |
25 |
27 #include <vector> |
26 #include <vector> |
|
27 #include <limits> |
|
28 |
|
29 #include <lemon/core.h> |
|
30 #include <lemon/maps.h> |
|
31 #include <lemon/path.h> |
|
32 #include <lemon/math.h> |
|
33 #include <lemon/static_graph.h> |
28 #include <lemon/adaptors.h> |
34 #include <lemon/adaptors.h> |
29 #include <lemon/path.h> |
|
30 |
|
31 #include <lemon/circulation.h> |
35 #include <lemon/circulation.h> |
32 #include <lemon/bellman_ford.h> |
36 #include <lemon/bellman_ford.h> |
33 #include <lemon/howard.h> |
37 #include <lemon/howard.h> |
34 |
38 |
35 namespace lemon { |
39 namespace lemon { |
36 |
40 |
37 /// \addtogroup min_cost_flow |
41 /// \addtogroup min_cost_flow_algs |
38 /// @{ |
42 /// @{ |
39 |
43 |
40 /// \brief Implementation of a cycle-canceling algorithm for |
44 /// \brief Implementation of cycle-canceling algorithms for |
41 /// finding a minimum cost flow. |
45 /// finding a \ref min_cost_flow "minimum cost flow". |
42 /// |
46 /// |
43 /// \ref CycleCanceling implements a cycle-canceling algorithm for |
47 /// \ref CycleCanceling implements three different cycle-canceling |
44 /// finding a minimum cost flow. |
48 /// algorithms for finding a \ref min_cost_flow "minimum cost flow". |
|
49 /// The most efficent one (both theoretically and practically) |
|
50 /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm, |
|
51 /// thus it is the default method. |
|
52 /// It is strongly polynomial, but in practice, it is typically much |
|
53 /// slower than the scaling algorithms and NetworkSimplex. |
45 /// |
54 /// |
46 /// \tparam Digraph The digraph type the algorithm runs on. |
55 /// Most of the parameters of the problem (except for the digraph) |
47 /// \tparam LowerMap The type of the lower bound map. |
56 /// can be given using separate functions, and the algorithm can be |
48 /// \tparam CapacityMap The type of the capacity (upper bound) map. |
57 /// executed using the \ref run() function. If some parameters are not |
49 /// \tparam CostMap The type of the cost (length) map. |
58 /// specified, then default values will be used. |
50 /// \tparam SupplyMap The type of the supply map. |
|
51 /// |
59 /// |
52 /// \warning |
60 /// \tparam GR The digraph type the algorithm runs on. |
53 /// - Arc capacities and costs should be \e non-negative \e integers. |
61 /// \tparam V The number type used for flow amounts, capacity bounds |
54 /// - Supply values should be \e signed \e integers. |
62 /// and supply values in the algorithm. By default, it is \c int. |
55 /// - The value types of the maps should be convertible to each other. |
63 /// \tparam C The number type used for costs and potentials in the |
56 /// - \c CostMap::Value must be signed type. |
64 /// algorithm. By default, it is the same as \c V. |
57 /// |
65 /// |
58 /// \note By default the \ref BellmanFord "Bellman-Ford" algorithm is |
66 /// \warning Both number types must be signed and all input data must |
59 /// used for negative cycle detection with limited iteration number. |
67 /// be integer. |
60 /// However \ref CycleCanceling also provides the "Minimum Mean |
68 /// \warning This algorithm does not support negative costs for such |
61 /// Cycle-Canceling" algorithm, which is \e strongly \e polynomial, |
69 /// arcs that have infinite upper bound. |
62 /// but rather slower in practice. |
|
63 /// To use this version of the algorithm, call \ref run() with \c true |
|
64 /// parameter. |
|
65 /// |
70 /// |
66 /// \author Peter Kovacs |
71 /// \note For more information about the three available methods, |
67 template < typename Digraph, |
72 /// see \ref Method. |
68 typename LowerMap = typename Digraph::template ArcMap<int>, |
73 #ifdef DOXYGEN |
69 typename CapacityMap = typename Digraph::template ArcMap<int>, |
74 template <typename GR, typename V, typename C> |
70 typename CostMap = typename Digraph::template ArcMap<int>, |
75 #else |
71 typename SupplyMap = typename Digraph::template NodeMap<int> > |
76 template <typename GR, typename V = int, typename C = V> |
|
77 #endif |
72 class CycleCanceling |
78 class CycleCanceling |
73 { |
79 { |
74 TEMPLATE_DIGRAPH_TYPEDEFS(Digraph); |
|
75 |
|
76 typedef typename CapacityMap::Value Capacity; |
|
77 typedef typename CostMap::Value Cost; |
|
78 typedef typename SupplyMap::Value Supply; |
|
79 typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap; |
|
80 typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap; |
|
81 |
|
82 typedef ResidualDigraph< const Digraph, |
|
83 CapacityArcMap, CapacityArcMap > ResDigraph; |
|
84 typedef typename ResDigraph::Node ResNode; |
|
85 typedef typename ResDigraph::NodeIt ResNodeIt; |
|
86 typedef typename ResDigraph::Arc ResArc; |
|
87 typedef typename ResDigraph::ArcIt ResArcIt; |
|
88 |
|
89 public: |
80 public: |
90 |
81 |
91 /// The type of the flow map. |
82 /// The type of the digraph |
92 typedef typename Digraph::template ArcMap<Capacity> FlowMap; |
83 typedef GR Digraph; |
93 /// The type of the potential map. |
84 /// The type of the flow amounts, capacity bounds and supply values |
94 typedef typename Digraph::template NodeMap<Cost> PotentialMap; |
85 typedef V Value; |
|
86 /// The type of the arc costs |
|
87 typedef C Cost; |
|
88 |
|
89 public: |
|
90 |
|
91 /// \brief Problem type constants for the \c run() function. |
|
92 /// |
|
93 /// Enum type containing the problem type constants that can be |
|
94 /// returned by the \ref run() function of the algorithm. |
|
95 enum ProblemType { |
|
96 /// The problem has no feasible solution (flow). |
|
97 INFEASIBLE, |
|
98 /// The problem has optimal solution (i.e. it is feasible and |
|
99 /// bounded), and the algorithm has found optimal flow and node |
|
100 /// potentials (primal and dual solutions). |
|
101 OPTIMAL, |
|
102 /// The digraph contains an arc of negative cost and infinite |
|
103 /// upper bound. It means that the objective function is unbounded |
|
104 /// on that arc, however, note that it could actually be bounded |
|
105 /// over the feasible flows, but this algroithm cannot handle |
|
106 /// these cases. |
|
107 UNBOUNDED |
|
108 }; |
|
109 |
|
110 /// \brief Constants for selecting the used method. |
|
111 /// |
|
112 /// Enum type containing constants for selecting the used method |
|
113 /// for the \ref run() function. |
|
114 /// |
|
115 /// \ref CycleCanceling provides three different cycle-canceling |
|
116 /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" |
|
117 /// is used, which proved to be the most efficient and the most robust |
|
118 /// on various test inputs. |
|
119 /// However, the other methods can be selected using the \ref run() |
|
120 /// function with the proper parameter. |
|
121 enum Method { |
|
122 /// A simple cycle-canceling method, which uses the |
|
123 /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration |
|
124 /// number for detecting negative cycles in the residual network. |
|
125 SIMPLE_CYCLE_CANCELING, |
|
126 /// The "Minimum Mean Cycle-Canceling" algorithm, which is a |
|
127 /// well-known strongly polynomial method. It improves along a |
|
128 /// \ref min_mean_cycle "minimum mean cycle" in each iteration. |
|
129 /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)). |
|
130 MINIMUM_MEAN_CYCLE_CANCELING, |
|
131 /// The "Cancel And Tighten" algorithm, which can be viewed as an |
|
132 /// improved version of the previous method. |
|
133 /// It is faster both in theory and in practice, its running time |
|
134 /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)). |
|
135 CANCEL_AND_TIGHTEN |
|
136 }; |
95 |
137 |
96 private: |
138 private: |
97 |
139 |
98 /// \brief Map adaptor class for handling residual arc costs. |
140 TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
99 /// |
141 |
100 /// Map adaptor class for handling residual arc costs. |
142 typedef std::vector<int> IntVector; |
101 class ResidualCostMap : public MapBase<ResArc, Cost> |
143 typedef std::vector<char> CharVector; |
|
144 typedef std::vector<double> DoubleVector; |
|
145 typedef std::vector<Value> ValueVector; |
|
146 typedef std::vector<Cost> CostVector; |
|
147 |
|
148 private: |
|
149 |
|
150 template <typename KT, typename VT> |
|
151 class VectorMap { |
|
152 public: |
|
153 typedef KT Key; |
|
154 typedef VT Value; |
|
155 |
|
156 VectorMap(std::vector<Value>& v) : _v(v) {} |
|
157 |
|
158 const Value& operator[](const Key& key) const { |
|
159 return _v[StaticDigraph::id(key)]; |
|
160 } |
|
161 |
|
162 Value& operator[](const Key& key) { |
|
163 return _v[StaticDigraph::id(key)]; |
|
164 } |
|
165 |
|
166 void set(const Key& key, const Value& val) { |
|
167 _v[StaticDigraph::id(key)] = val; |
|
168 } |
|
169 |
|
170 private: |
|
171 std::vector<Value>& _v; |
|
172 }; |
|
173 |
|
174 typedef VectorMap<StaticDigraph::Node, Cost> CostNodeMap; |
|
175 typedef VectorMap<StaticDigraph::Arc, Cost> CostArcMap; |
|
176 |
|
177 private: |
|
178 |
|
179 |
|
180 // Data related to the underlying digraph |
|
181 const GR &_graph; |
|
182 int _node_num; |
|
183 int _arc_num; |
|
184 int _res_node_num; |
|
185 int _res_arc_num; |
|
186 int _root; |
|
187 |
|
188 // Parameters of the problem |
|
189 bool _have_lower; |
|
190 Value _sum_supply; |
|
191 |
|
192 // Data structures for storing the digraph |
|
193 IntNodeMap _node_id; |
|
194 IntArcMap _arc_idf; |
|
195 IntArcMap _arc_idb; |
|
196 IntVector _first_out; |
|
197 CharVector _forward; |
|
198 IntVector _source; |
|
199 IntVector _target; |
|
200 IntVector _reverse; |
|
201 |
|
202 // Node and arc data |
|
203 ValueVector _lower; |
|
204 ValueVector _upper; |
|
205 CostVector _cost; |
|
206 ValueVector _supply; |
|
207 |
|
208 ValueVector _res_cap; |
|
209 CostVector _pi; |
|
210 |
|
211 // Data for a StaticDigraph structure |
|
212 typedef std::pair<int, int> IntPair; |
|
213 StaticDigraph _sgr; |
|
214 std::vector<IntPair> _arc_vec; |
|
215 std::vector<Cost> _cost_vec; |
|
216 IntVector _id_vec; |
|
217 CostArcMap _cost_map; |
|
218 CostNodeMap _pi_map; |
|
219 |
|
220 public: |
|
221 |
|
222 /// \brief Constant for infinite upper bounds (capacities). |
|
223 /// |
|
224 /// Constant for infinite upper bounds (capacities). |
|
225 /// It is \c std::numeric_limits<Value>::infinity() if available, |
|
226 /// \c std::numeric_limits<Value>::max() otherwise. |
|
227 const Value INF; |
|
228 |
|
229 public: |
|
230 |
|
231 /// \brief Constructor. |
|
232 /// |
|
233 /// The constructor of the class. |
|
234 /// |
|
235 /// \param graph The digraph the algorithm runs on. |
|
236 CycleCanceling(const GR& graph) : |
|
237 _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph), |
|
238 _cost_map(_cost_vec), _pi_map(_pi), |
|
239 INF(std::numeric_limits<Value>::has_infinity ? |
|
240 std::numeric_limits<Value>::infinity() : |
|
241 std::numeric_limits<Value>::max()) |
102 { |
242 { |
103 private: |
243 // Check the number types |
104 |
244 LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
105 const CostMap &_cost_map; |
245 "The flow type of CycleCanceling must be signed"); |
106 |
246 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
107 public: |
247 "The cost type of CycleCanceling must be signed"); |
108 |
248 |
109 ///\e |
249 // Resize vectors |
110 ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {} |
250 _node_num = countNodes(_graph); |
111 |
251 _arc_num = countArcs(_graph); |
112 ///\e |
252 _res_node_num = _node_num + 1; |
113 Cost operator[](const ResArc &e) const { |
253 _res_arc_num = 2 * (_arc_num + _node_num); |
114 return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e]; |
254 _root = _node_num; |
115 } |
255 |
116 |
256 _first_out.resize(_res_node_num + 1); |
117 }; //class ResidualCostMap |
257 _forward.resize(_res_arc_num); |
118 |
258 _source.resize(_res_arc_num); |
119 private: |
259 _target.resize(_res_arc_num); |
120 |
260 _reverse.resize(_res_arc_num); |
121 // The maximum number of iterations for the first execution of the |
261 |
122 // Bellman-Ford algorithm. It should be at least 2. |
262 _lower.resize(_res_arc_num); |
123 static const int BF_FIRST_LIMIT = 2; |
263 _upper.resize(_res_arc_num); |
124 // The iteration limit for the Bellman-Ford algorithm is multiplied |
264 _cost.resize(_res_arc_num); |
125 // by BF_LIMIT_FACTOR/100 in every round. |
265 _supply.resize(_res_node_num); |
126 static const int BF_LIMIT_FACTOR = 150; |
266 |
127 |
267 _res_cap.resize(_res_arc_num); |
128 private: |
268 _pi.resize(_res_node_num); |
129 |
269 |
130 // The digraph the algorithm runs on |
270 _arc_vec.reserve(_res_arc_num); |
131 const Digraph &_graph; |
271 _cost_vec.reserve(_res_arc_num); |
132 // The original lower bound map |
272 _id_vec.reserve(_res_arc_num); |
133 const LowerMap *_lower; |
273 |
134 // The modified capacity map |
274 // Copy the graph |
135 CapacityArcMap _capacity; |
275 int i = 0, j = 0, k = 2 * _arc_num + _node_num; |
136 // The original cost map |
276 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
137 const CostMap &_cost; |
277 _node_id[n] = i; |
138 // The modified supply map |
278 } |
139 SupplyNodeMap _supply; |
279 i = 0; |
140 bool _valid_supply; |
280 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
141 |
281 _first_out[i] = j; |
142 // Arc map of the current flow |
282 for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
143 FlowMap *_flow; |
283 _arc_idf[a] = j; |
144 bool _local_flow; |
284 _forward[j] = true; |
145 // Node map of the current potentials |
285 _source[j] = i; |
146 PotentialMap *_potential; |
286 _target[j] = _node_id[_graph.runningNode(a)]; |
147 bool _local_potential; |
287 } |
148 |
288 for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) { |
149 // The residual digraph |
289 _arc_idb[a] = j; |
150 ResDigraph *_res_graph; |
290 _forward[j] = false; |
151 // The residual cost map |
291 _source[j] = i; |
152 ResidualCostMap _res_cost; |
292 _target[j] = _node_id[_graph.runningNode(a)]; |
153 |
293 } |
154 public: |
294 _forward[j] = false; |
155 |
295 _source[j] = i; |
156 /// \brief General constructor (with lower bounds). |
296 _target[j] = _root; |
157 /// |
297 _reverse[j] = k; |
158 /// General constructor (with lower bounds). |
298 _forward[k] = true; |
159 /// |
299 _source[k] = _root; |
160 /// \param digraph The digraph the algorithm runs on. |
300 _target[k] = i; |
161 /// \param lower The lower bounds of the arcs. |
301 _reverse[k] = j; |
162 /// \param capacity The capacities (upper bounds) of the arcs. |
302 ++j; ++k; |
163 /// \param cost The cost (length) values of the arcs. |
303 } |
164 /// \param supply The supply values of the nodes (signed). |
304 _first_out[i] = j; |
165 CycleCanceling( const Digraph &digraph, |
305 _first_out[_res_node_num] = k; |
166 const LowerMap &lower, |
306 for (ArcIt a(_graph); a != INVALID; ++a) { |
167 const CapacityMap &capacity, |
307 int fi = _arc_idf[a]; |
168 const CostMap &cost, |
308 int bi = _arc_idb[a]; |
169 const SupplyMap &supply ) : |
309 _reverse[fi] = bi; |
170 _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost), |
310 _reverse[bi] = fi; |
171 _supply(digraph), _flow(NULL), _local_flow(false), |
311 } |
172 _potential(NULL), _local_potential(false), |
312 |
173 _res_graph(NULL), _res_cost(_cost) |
313 // Reset parameters |
174 { |
314 reset(); |
175 // Check the sum of supply values |
315 } |
176 Supply sum = 0; |
316 |
|
317 /// \name Parameters |
|
318 /// The parameters of the algorithm can be specified using these |
|
319 /// functions. |
|
320 |
|
321 /// @{ |
|
322 |
|
323 /// \brief Set the lower bounds on the arcs. |
|
324 /// |
|
325 /// This function sets the lower bounds on the arcs. |
|
326 /// If it is not used before calling \ref run(), the lower bounds |
|
327 /// will be set to zero on all arcs. |
|
328 /// |
|
329 /// \param map An arc map storing the lower bounds. |
|
330 /// Its \c Value type must be convertible to the \c Value type |
|
331 /// of the algorithm. |
|
332 /// |
|
333 /// \return <tt>(*this)</tt> |
|
334 template <typename LowerMap> |
|
335 CycleCanceling& lowerMap(const LowerMap& map) { |
|
336 _have_lower = true; |
|
337 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
338 _lower[_arc_idf[a]] = map[a]; |
|
339 _lower[_arc_idb[a]] = map[a]; |
|
340 } |
|
341 return *this; |
|
342 } |
|
343 |
|
344 /// \brief Set the upper bounds (capacities) on the arcs. |
|
345 /// |
|
346 /// This function sets the upper bounds (capacities) on the arcs. |
|
347 /// If it is not used before calling \ref run(), the upper bounds |
|
348 /// will be set to \ref INF on all arcs (i.e. the flow value will be |
|
349 /// unbounded from above). |
|
350 /// |
|
351 /// \param map An arc map storing the upper bounds. |
|
352 /// Its \c Value type must be convertible to the \c Value type |
|
353 /// of the algorithm. |
|
354 /// |
|
355 /// \return <tt>(*this)</tt> |
|
356 template<typename UpperMap> |
|
357 CycleCanceling& upperMap(const UpperMap& map) { |
|
358 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
359 _upper[_arc_idf[a]] = map[a]; |
|
360 } |
|
361 return *this; |
|
362 } |
|
363 |
|
364 /// \brief Set the costs of the arcs. |
|
365 /// |
|
366 /// This function sets the costs of the arcs. |
|
367 /// If it is not used before calling \ref run(), the costs |
|
368 /// will be set to \c 1 on all arcs. |
|
369 /// |
|
370 /// \param map An arc map storing the costs. |
|
371 /// Its \c Value type must be convertible to the \c Cost type |
|
372 /// of the algorithm. |
|
373 /// |
|
374 /// \return <tt>(*this)</tt> |
|
375 template<typename CostMap> |
|
376 CycleCanceling& costMap(const CostMap& map) { |
|
377 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
378 _cost[_arc_idf[a]] = map[a]; |
|
379 _cost[_arc_idb[a]] = -map[a]; |
|
380 } |
|
381 return *this; |
|
382 } |
|
383 |
|
384 /// \brief Set the supply values of the nodes. |
|
385 /// |
|
386 /// This function sets the supply values of the nodes. |
|
387 /// If neither this function nor \ref stSupply() is used before |
|
388 /// calling \ref run(), the supply of each node will be set to zero. |
|
389 /// |
|
390 /// \param map A node map storing the supply values. |
|
391 /// Its \c Value type must be convertible to the \c Value type |
|
392 /// of the algorithm. |
|
393 /// |
|
394 /// \return <tt>(*this)</tt> |
|
395 template<typename SupplyMap> |
|
396 CycleCanceling& supplyMap(const SupplyMap& map) { |
177 for (NodeIt n(_graph); n != INVALID; ++n) { |
397 for (NodeIt n(_graph); n != INVALID; ++n) { |
178 _supply[n] = supply[n]; |
398 _supply[_node_id[n]] = map[n]; |
179 sum += _supply[n]; |
399 } |
180 } |
400 return *this; |
181 _valid_supply = sum == 0; |
401 } |
182 |
402 |
183 // Remove non-zero lower bounds |
403 /// \brief Set single source and target nodes and a supply value. |
184 for (ArcIt e(_graph); e != INVALID; ++e) { |
404 /// |
185 _capacity[e] = capacity[e]; |
405 /// This function sets a single source node and a single target node |
186 if (lower[e] != 0) { |
406 /// and the required flow value. |
187 _capacity[e] -= lower[e]; |
407 /// If neither this function nor \ref supplyMap() is used before |
188 _supply[_graph.source(e)] -= lower[e]; |
408 /// calling \ref run(), the supply of each node will be set to zero. |
189 _supply[_graph.target(e)] += lower[e]; |
409 /// |
190 } |
410 /// Using this function has the same effect as using \ref supplyMap() |
191 } |
411 /// with such a map in which \c k is assigned to \c s, \c -k is |
192 } |
412 /// assigned to \c t and all other nodes have zero supply value. |
193 /* |
413 /// |
194 /// \brief General constructor (without lower bounds). |
|
195 /// |
|
196 /// General constructor (without lower bounds). |
|
197 /// |
|
198 /// \param digraph The digraph the algorithm runs on. |
|
199 /// \param capacity The capacities (upper bounds) of the arcs. |
|
200 /// \param cost The cost (length) values of the arcs. |
|
201 /// \param supply The supply values of the nodes (signed). |
|
202 CycleCanceling( const Digraph &digraph, |
|
203 const CapacityMap &capacity, |
|
204 const CostMap &cost, |
|
205 const SupplyMap &supply ) : |
|
206 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost), |
|
207 _supply(supply), _flow(NULL), _local_flow(false), |
|
208 _potential(NULL), _local_potential(false), _res_graph(NULL), |
|
209 _res_cost(_cost) |
|
210 { |
|
211 // Check the sum of supply values |
|
212 Supply sum = 0; |
|
213 for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n]; |
|
214 _valid_supply = sum == 0; |
|
215 } |
|
216 |
|
217 /// \brief Simple constructor (with lower bounds). |
|
218 /// |
|
219 /// Simple constructor (with lower bounds). |
|
220 /// |
|
221 /// \param digraph The digraph the algorithm runs on. |
|
222 /// \param lower The lower bounds of the arcs. |
|
223 /// \param capacity The capacities (upper bounds) of the arcs. |
|
224 /// \param cost The cost (length) values of the arcs. |
|
225 /// \param s The source node. |
414 /// \param s The source node. |
226 /// \param t The target node. |
415 /// \param t The target node. |
227 /// \param flow_value The required amount of flow from node \c s |
416 /// \param k The required amount of flow from node \c s to node \c t |
228 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
417 /// (i.e. the supply of \c s and the demand of \c t). |
229 CycleCanceling( const Digraph &digraph, |
418 /// |
230 const LowerMap &lower, |
419 /// \return <tt>(*this)</tt> |
231 const CapacityMap &capacity, |
420 CycleCanceling& stSupply(const Node& s, const Node& t, Value k) { |
232 const CostMap &cost, |
421 for (int i = 0; i != _res_node_num; ++i) { |
233 Node s, Node t, |
422 _supply[i] = 0; |
234 Supply flow_value ) : |
423 } |
235 _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost), |
424 _supply[_node_id[s]] = k; |
236 _supply(digraph, 0), _flow(NULL), _local_flow(false), |
425 _supply[_node_id[t]] = -k; |
237 _potential(NULL), _local_potential(false), _res_graph(NULL), |
|
238 _res_cost(_cost) |
|
239 { |
|
240 // Remove non-zero lower bounds |
|
241 _supply[s] = flow_value; |
|
242 _supply[t] = -flow_value; |
|
243 for (ArcIt e(_graph); e != INVALID; ++e) { |
|
244 if (lower[e] != 0) { |
|
245 _capacity[e] -= lower[e]; |
|
246 _supply[_graph.source(e)] -= lower[e]; |
|
247 _supply[_graph.target(e)] += lower[e]; |
|
248 } |
|
249 } |
|
250 _valid_supply = true; |
|
251 } |
|
252 |
|
253 /// \brief Simple constructor (without lower bounds). |
|
254 /// |
|
255 /// Simple constructor (without lower bounds). |
|
256 /// |
|
257 /// \param digraph The digraph the algorithm runs on. |
|
258 /// \param capacity The capacities (upper bounds) of the arcs. |
|
259 /// \param cost The cost (length) values of the arcs. |
|
260 /// \param s The source node. |
|
261 /// \param t The target node. |
|
262 /// \param flow_value The required amount of flow from node \c s |
|
263 /// to node \c t (i.e. the supply of \c s and the demand of \c t). |
|
264 CycleCanceling( const Digraph &digraph, |
|
265 const CapacityMap &capacity, |
|
266 const CostMap &cost, |
|
267 Node s, Node t, |
|
268 Supply flow_value ) : |
|
269 _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost), |
|
270 _supply(digraph, 0), _flow(NULL), _local_flow(false), |
|
271 _potential(NULL), _local_potential(false), _res_graph(NULL), |
|
272 _res_cost(_cost) |
|
273 { |
|
274 _supply[s] = flow_value; |
|
275 _supply[t] = -flow_value; |
|
276 _valid_supply = true; |
|
277 } |
|
278 */ |
|
279 /// Destructor. |
|
280 ~CycleCanceling() { |
|
281 if (_local_flow) delete _flow; |
|
282 if (_local_potential) delete _potential; |
|
283 delete _res_graph; |
|
284 } |
|
285 |
|
286 /// \brief Set the flow map. |
|
287 /// |
|
288 /// Set the flow map. |
|
289 /// |
|
290 /// \return \c (*this) |
|
291 CycleCanceling& flowMap(FlowMap &map) { |
|
292 if (_local_flow) { |
|
293 delete _flow; |
|
294 _local_flow = false; |
|
295 } |
|
296 _flow = ↦ |
|
297 return *this; |
426 return *this; |
298 } |
427 } |
299 |
428 |
300 /// \brief Set the potential map. |
429 /// @} |
301 /// |
430 |
302 /// Set the potential map. |
431 /// \name Execution control |
303 /// |
432 /// The algorithm can be executed using \ref run(). |
304 /// \return \c (*this) |
433 |
305 CycleCanceling& potentialMap(PotentialMap &map) { |
434 /// @{ |
306 if (_local_potential) { |
435 |
307 delete _potential; |
436 /// \brief Run the algorithm. |
308 _local_potential = false; |
437 /// |
309 } |
438 /// This function runs the algorithm. |
310 _potential = ↦ |
439 /// The paramters can be specified using functions \ref lowerMap(), |
|
440 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(). |
|
441 /// For example, |
|
442 /// \code |
|
443 /// CycleCanceling<ListDigraph> cc(graph); |
|
444 /// cc.lowerMap(lower).upperMap(upper).costMap(cost) |
|
445 /// .supplyMap(sup).run(); |
|
446 /// \endcode |
|
447 /// |
|
448 /// This function can be called more than once. All the parameters |
|
449 /// that have been given are kept for the next call, unless |
|
450 /// \ref reset() is called, thus only the modified parameters |
|
451 /// have to be set again. See \ref reset() for examples. |
|
452 /// However, the underlying digraph must not be modified after this |
|
453 /// class have been constructed, since it copies and extends the graph. |
|
454 /// |
|
455 /// \param method The cycle-canceling method that will be used. |
|
456 /// For more information, see \ref Method. |
|
457 /// |
|
458 /// \return \c INFEASIBLE if no feasible flow exists, |
|
459 /// \n \c OPTIMAL if the problem has optimal solution |
|
460 /// (i.e. it is feasible and bounded), and the algorithm has found |
|
461 /// optimal flow and node potentials (primal and dual solutions), |
|
462 /// \n \c UNBOUNDED if the digraph contains an arc of negative cost |
|
463 /// and infinite upper bound. It means that the objective function |
|
464 /// is unbounded on that arc, however, note that it could actually be |
|
465 /// bounded over the feasible flows, but this algroithm cannot handle |
|
466 /// these cases. |
|
467 /// |
|
468 /// \see ProblemType, Method |
|
469 ProblemType run(Method method = CANCEL_AND_TIGHTEN) { |
|
470 ProblemType pt = init(); |
|
471 if (pt != OPTIMAL) return pt; |
|
472 start(method); |
|
473 return OPTIMAL; |
|
474 } |
|
475 |
|
476 /// \brief Reset all the parameters that have been given before. |
|
477 /// |
|
478 /// This function resets all the paramaters that have been given |
|
479 /// before using functions \ref lowerMap(), \ref upperMap(), |
|
480 /// \ref costMap(), \ref supplyMap(), \ref stSupply(). |
|
481 /// |
|
482 /// It is useful for multiple run() calls. If this function is not |
|
483 /// used, all the parameters given before are kept for the next |
|
484 /// \ref run() call. |
|
485 /// However, the underlying digraph must not be modified after this |
|
486 /// class have been constructed, since it copies and extends the graph. |
|
487 /// |
|
488 /// For example, |
|
489 /// \code |
|
490 /// CycleCanceling<ListDigraph> cs(graph); |
|
491 /// |
|
492 /// // First run |
|
493 /// cc.lowerMap(lower).upperMap(upper).costMap(cost) |
|
494 /// .supplyMap(sup).run(); |
|
495 /// |
|
496 /// // Run again with modified cost map (reset() is not called, |
|
497 /// // so only the cost map have to be set again) |
|
498 /// cost[e] += 100; |
|
499 /// cc.costMap(cost).run(); |
|
500 /// |
|
501 /// // Run again from scratch using reset() |
|
502 /// // (the lower bounds will be set to zero on all arcs) |
|
503 /// cc.reset(); |
|
504 /// cc.upperMap(capacity).costMap(cost) |
|
505 /// .supplyMap(sup).run(); |
|
506 /// \endcode |
|
507 /// |
|
508 /// \return <tt>(*this)</tt> |
|
509 CycleCanceling& reset() { |
|
510 for (int i = 0; i != _res_node_num; ++i) { |
|
511 _supply[i] = 0; |
|
512 } |
|
513 int limit = _first_out[_root]; |
|
514 for (int j = 0; j != limit; ++j) { |
|
515 _lower[j] = 0; |
|
516 _upper[j] = INF; |
|
517 _cost[j] = _forward[j] ? 1 : -1; |
|
518 } |
|
519 for (int j = limit; j != _res_arc_num; ++j) { |
|
520 _lower[j] = 0; |
|
521 _upper[j] = INF; |
|
522 _cost[j] = 0; |
|
523 _cost[_reverse[j]] = 0; |
|
524 } |
|
525 _have_lower = false; |
311 return *this; |
526 return *this; |
312 } |
527 } |
313 |
528 |
314 /// \name Execution control |
529 /// @} |
|
530 |
|
531 /// \name Query Functions |
|
532 /// The results of the algorithm can be obtained using these |
|
533 /// functions.\n |
|
534 /// The \ref run() function must be called before using them. |
315 |
535 |
316 /// @{ |
536 /// @{ |
317 |
537 |
318 /// \brief Run the algorithm. |
538 /// \brief Return the total cost of the found flow. |
319 /// |
539 /// |
320 /// Run the algorithm. |
540 /// This function returns the total cost of the found flow. |
321 /// |
541 /// Its complexity is O(e). |
322 /// \param min_mean_cc Set this parameter to \c true to run the |
542 /// |
323 /// "Minimum Mean Cycle-Canceling" algorithm, which is strongly |
543 /// \note The return type of the function can be specified as a |
324 /// polynomial, but rather slower in practice. |
544 /// template parameter. For example, |
325 /// |
545 /// \code |
326 /// \return \c true if a feasible flow can be found. |
546 /// cc.totalCost<double>(); |
327 bool run(bool min_mean_cc = false) { |
547 /// \endcode |
328 return init() && start(min_mean_cc); |
548 /// It is useful if the total cost cannot be stored in the \c Cost |
|
549 /// type of the algorithm, which is the default return type of the |
|
550 /// function. |
|
551 /// |
|
552 /// \pre \ref run() must be called before using this function. |
|
553 template <typename Number> |
|
554 Number totalCost() const { |
|
555 Number c = 0; |
|
556 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
557 int i = _arc_idb[a]; |
|
558 c += static_cast<Number>(_res_cap[i]) * |
|
559 (-static_cast<Number>(_cost[i])); |
|
560 } |
|
561 return c; |
|
562 } |
|
563 |
|
564 #ifndef DOXYGEN |
|
565 Cost totalCost() const { |
|
566 return totalCost<Cost>(); |
|
567 } |
|
568 #endif |
|
569 |
|
570 /// \brief Return the flow on the given arc. |
|
571 /// |
|
572 /// This function returns the flow on the given arc. |
|
573 /// |
|
574 /// \pre \ref run() must be called before using this function. |
|
575 Value flow(const Arc& a) const { |
|
576 return _res_cap[_arc_idb[a]]; |
|
577 } |
|
578 |
|
579 /// \brief Return the flow map (the primal solution). |
|
580 /// |
|
581 /// This function copies the flow value on each arc into the given |
|
582 /// map. The \c Value type of the algorithm must be convertible to |
|
583 /// the \c Value type of the map. |
|
584 /// |
|
585 /// \pre \ref run() must be called before using this function. |
|
586 template <typename FlowMap> |
|
587 void flowMap(FlowMap &map) const { |
|
588 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
589 map.set(a, _res_cap[_arc_idb[a]]); |
|
590 } |
|
591 } |
|
592 |
|
593 /// \brief Return the potential (dual value) of the given node. |
|
594 /// |
|
595 /// This function returns the potential (dual value) of the |
|
596 /// given node. |
|
597 /// |
|
598 /// \pre \ref run() must be called before using this function. |
|
599 Cost potential(const Node& n) const { |
|
600 return static_cast<Cost>(_pi[_node_id[n]]); |
|
601 } |
|
602 |
|
603 /// \brief Return the potential map (the dual solution). |
|
604 /// |
|
605 /// This function copies the potential (dual value) of each node |
|
606 /// into the given map. |
|
607 /// The \c Cost type of the algorithm must be convertible to the |
|
608 /// \c Value type of the map. |
|
609 /// |
|
610 /// \pre \ref run() must be called before using this function. |
|
611 template <typename PotentialMap> |
|
612 void potentialMap(PotentialMap &map) const { |
|
613 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
614 map.set(n, static_cast<Cost>(_pi[_node_id[n]])); |
|
615 } |
329 } |
616 } |
330 |
617 |
331 /// @} |
618 /// @} |
332 |
619 |
333 /// \name Query Functions |
|
334 /// The result of the algorithm can be obtained using these |
|
335 /// functions.\n |
|
336 /// \ref lemon::CycleCanceling::run() "run()" must be called before |
|
337 /// using them. |
|
338 |
|
339 /// @{ |
|
340 |
|
341 /// \brief Return a const reference to the arc map storing the |
|
342 /// found flow. |
|
343 /// |
|
344 /// Return a const reference to the arc map storing the found flow. |
|
345 /// |
|
346 /// \pre \ref run() must be called before using this function. |
|
347 const FlowMap& flowMap() const { |
|
348 return *_flow; |
|
349 } |
|
350 |
|
351 /// \brief Return a const reference to the node map storing the |
|
352 /// found potentials (the dual solution). |
|
353 /// |
|
354 /// Return a const reference to the node map storing the found |
|
355 /// potentials (the dual solution). |
|
356 /// |
|
357 /// \pre \ref run() must be called before using this function. |
|
358 const PotentialMap& potentialMap() const { |
|
359 return *_potential; |
|
360 } |
|
361 |
|
362 /// \brief Return the flow on the given arc. |
|
363 /// |
|
364 /// Return the flow on the given arc. |
|
365 /// |
|
366 /// \pre \ref run() must be called before using this function. |
|
367 Capacity flow(const Arc& arc) const { |
|
368 return (*_flow)[arc]; |
|
369 } |
|
370 |
|
371 /// \brief Return the potential of the given node. |
|
372 /// |
|
373 /// Return the potential of the given node. |
|
374 /// |
|
375 /// \pre \ref run() must be called before using this function. |
|
376 Cost potential(const Node& node) const { |
|
377 return (*_potential)[node]; |
|
378 } |
|
379 |
|
380 /// \brief Return the total cost of the found flow. |
|
381 /// |
|
382 /// Return the total cost of the found flow. The complexity of the |
|
383 /// function is \f$ O(e) \f$. |
|
384 /// |
|
385 /// \pre \ref run() must be called before using this function. |
|
386 Cost totalCost() const { |
|
387 Cost c = 0; |
|
388 for (ArcIt e(_graph); e != INVALID; ++e) |
|
389 c += (*_flow)[e] * _cost[e]; |
|
390 return c; |
|
391 } |
|
392 |
|
393 /// @} |
|
394 |
|
395 private: |
620 private: |
396 |
621 |
397 /// Initialize the algorithm. |
622 // Initialize the algorithm |
398 bool init() { |
623 ProblemType init() { |
399 if (!_valid_supply) return false; |
624 if (_res_node_num <= 1) return INFEASIBLE; |
400 |
625 |
401 // Initializing flow and potential maps |
626 // Check the sum of supply values |
402 if (!_flow) { |
627 _sum_supply = 0; |
403 _flow = new FlowMap(_graph); |
628 for (int i = 0; i != _root; ++i) { |
404 _local_flow = true; |
629 _sum_supply += _supply[i]; |
405 } |
630 } |
406 if (!_potential) { |
631 if (_sum_supply > 0) return INFEASIBLE; |
407 _potential = new PotentialMap(_graph); |
632 |
408 _local_potential = true; |
633 |
409 } |
634 // Initialize vectors |
410 |
635 for (int i = 0; i != _res_node_num; ++i) { |
411 _res_graph = new ResDigraph(_graph, _capacity, *_flow); |
636 _pi[i] = 0; |
412 |
637 } |
413 // Finding a feasible flow using Circulation |
638 ValueVector excess(_supply); |
414 Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap, |
639 |
415 SupplyMap > |
640 // Remove infinite upper bounds and check negative arcs |
416 circulation( _graph, constMap<Arc>(Capacity(0)), _capacity, |
641 const Value MAX = std::numeric_limits<Value>::max(); |
417 _supply ); |
642 int last_out; |
418 return circulation.flowMap(*_flow).run(); |
643 if (_have_lower) { |
419 } |
644 for (int i = 0; i != _root; ++i) { |
420 |
645 last_out = _first_out[i+1]; |
421 bool start(bool min_mean_cc) { |
646 for (int j = _first_out[i]; j != last_out; ++j) { |
422 if (min_mean_cc) |
647 if (_forward[j]) { |
423 startMinMean(); |
648 Value c = _cost[j] < 0 ? _upper[j] : _lower[j]; |
424 else |
649 if (c >= MAX) return UNBOUNDED; |
425 start(); |
650 excess[i] -= c; |
426 |
651 excess[_target[j]] += c; |
427 // Handling non-zero lower bounds |
652 } |
428 if (_lower) { |
653 } |
429 for (ArcIt e(_graph); e != INVALID; ++e) |
654 } |
430 (*_flow)[e] += (*_lower)[e]; |
655 } else { |
431 } |
656 for (int i = 0; i != _root; ++i) { |
432 return true; |
657 last_out = _first_out[i+1]; |
433 } |
658 for (int j = _first_out[i]; j != last_out; ++j) { |
434 |
659 if (_forward[j] && _cost[j] < 0) { |
435 /// \brief Execute the algorithm using \ref BellmanFord. |
660 Value c = _upper[j]; |
436 /// |
661 if (c >= MAX) return UNBOUNDED; |
437 /// Execute the algorithm using the \ref BellmanFord |
662 excess[i] -= c; |
438 /// "Bellman-Ford" algorithm for negative cycle detection with |
663 excess[_target[j]] += c; |
439 /// successively larger limit for the number of iterations. |
664 } |
440 void start() { |
665 } |
441 typename BellmanFord<ResDigraph, ResidualCostMap>::PredMap pred(*_res_graph); |
666 } |
442 typename ResDigraph::template NodeMap<int> visited(*_res_graph); |
667 } |
443 std::vector<ResArc> cycle; |
668 Value ex, max_cap = 0; |
444 int node_num = countNodes(_graph); |
669 for (int i = 0; i != _res_node_num; ++i) { |
|
670 ex = excess[i]; |
|
671 if (ex < 0) max_cap -= ex; |
|
672 } |
|
673 for (int j = 0; j != _res_arc_num; ++j) { |
|
674 if (_upper[j] >= MAX) _upper[j] = max_cap; |
|
675 } |
|
676 |
|
677 // Initialize maps for Circulation and remove non-zero lower bounds |
|
678 ConstMap<Arc, Value> low(0); |
|
679 typedef typename Digraph::template ArcMap<Value> ValueArcMap; |
|
680 typedef typename Digraph::template NodeMap<Value> ValueNodeMap; |
|
681 ValueArcMap cap(_graph), flow(_graph); |
|
682 ValueNodeMap sup(_graph); |
|
683 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
684 sup[n] = _supply[_node_id[n]]; |
|
685 } |
|
686 if (_have_lower) { |
|
687 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
688 int j = _arc_idf[a]; |
|
689 Value c = _lower[j]; |
|
690 cap[a] = _upper[j] - c; |
|
691 sup[_graph.source(a)] -= c; |
|
692 sup[_graph.target(a)] += c; |
|
693 } |
|
694 } else { |
|
695 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
696 cap[a] = _upper[_arc_idf[a]]; |
|
697 } |
|
698 } |
|
699 |
|
700 // Find a feasible flow using Circulation |
|
701 Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap> |
|
702 circ(_graph, low, cap, sup); |
|
703 if (!circ.flowMap(flow).run()) return INFEASIBLE; |
|
704 |
|
705 // Set residual capacities and handle GEQ supply type |
|
706 if (_sum_supply < 0) { |
|
707 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
708 Value fa = flow[a]; |
|
709 _res_cap[_arc_idf[a]] = cap[a] - fa; |
|
710 _res_cap[_arc_idb[a]] = fa; |
|
711 sup[_graph.source(a)] -= fa; |
|
712 sup[_graph.target(a)] += fa; |
|
713 } |
|
714 for (NodeIt n(_graph); n != INVALID; ++n) { |
|
715 excess[_node_id[n]] = sup[n]; |
|
716 } |
|
717 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
|
718 int u = _target[a]; |
|
719 int ra = _reverse[a]; |
|
720 _res_cap[a] = -_sum_supply + 1; |
|
721 _res_cap[ra] = -excess[u]; |
|
722 _cost[a] = 0; |
|
723 _cost[ra] = 0; |
|
724 } |
|
725 } else { |
|
726 for (ArcIt a(_graph); a != INVALID; ++a) { |
|
727 Value fa = flow[a]; |
|
728 _res_cap[_arc_idf[a]] = cap[a] - fa; |
|
729 _res_cap[_arc_idb[a]] = fa; |
|
730 } |
|
731 for (int a = _first_out[_root]; a != _res_arc_num; ++a) { |
|
732 int ra = _reverse[a]; |
|
733 _res_cap[a] = 1; |
|
734 _res_cap[ra] = 0; |
|
735 _cost[a] = 0; |
|
736 _cost[ra] = 0; |
|
737 } |
|
738 } |
|
739 |
|
740 return OPTIMAL; |
|
741 } |
|
742 |
|
743 // Build a StaticDigraph structure containing the current |
|
744 // residual network |
|
745 void buildResidualNetwork() { |
|
746 _arc_vec.clear(); |
|
747 _cost_vec.clear(); |
|
748 _id_vec.clear(); |
|
749 for (int j = 0; j != _res_arc_num; ++j) { |
|
750 if (_res_cap[j] > 0) { |
|
751 _arc_vec.push_back(IntPair(_source[j], _target[j])); |
|
752 _cost_vec.push_back(_cost[j]); |
|
753 _id_vec.push_back(j); |
|
754 } |
|
755 } |
|
756 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); |
|
757 } |
|
758 |
|
759 // Execute the algorithm and transform the results |
|
760 void start(Method method) { |
|
761 // Execute the algorithm |
|
762 switch (method) { |
|
763 case SIMPLE_CYCLE_CANCELING: |
|
764 startSimpleCycleCanceling(); |
|
765 break; |
|
766 case MINIMUM_MEAN_CYCLE_CANCELING: |
|
767 startMinMeanCycleCanceling(); |
|
768 break; |
|
769 case CANCEL_AND_TIGHTEN: |
|
770 startCancelAndTighten(); |
|
771 break; |
|
772 } |
|
773 |
|
774 // Compute node potentials |
|
775 if (method != SIMPLE_CYCLE_CANCELING) { |
|
776 buildResidualNetwork(); |
|
777 typename BellmanFord<StaticDigraph, CostArcMap> |
|
778 ::template SetDistMap<CostNodeMap>::Create bf(_sgr, _cost_map); |
|
779 bf.distMap(_pi_map); |
|
780 bf.init(0); |
|
781 bf.start(); |
|
782 } |
|
783 |
|
784 // Handle non-zero lower bounds |
|
785 if (_have_lower) { |
|
786 int limit = _first_out[_root]; |
|
787 for (int j = 0; j != limit; ++j) { |
|
788 if (!_forward[j]) _res_cap[j] += _lower[j]; |
|
789 } |
|
790 } |
|
791 } |
|
792 |
|
793 // Execute the "Simple Cycle Canceling" method |
|
794 void startSimpleCycleCanceling() { |
|
795 // Constants for computing the iteration limits |
|
796 const int BF_FIRST_LIMIT = 2; |
|
797 const double BF_LIMIT_FACTOR = 1.5; |
|
798 |
|
799 typedef VectorMap<StaticDigraph::Arc, Value> FilterMap; |
|
800 typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph; |
|
801 typedef VectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap; |
|
802 typedef typename BellmanFord<ResDigraph, CostArcMap> |
|
803 ::template SetDistMap<CostNodeMap> |
|
804 ::template SetPredMap<PredMap>::Create BF; |
|
805 |
|
806 // Build the residual network |
|
807 _arc_vec.clear(); |
|
808 _cost_vec.clear(); |
|
809 for (int j = 0; j != _res_arc_num; ++j) { |
|
810 _arc_vec.push_back(IntPair(_source[j], _target[j])); |
|
811 _cost_vec.push_back(_cost[j]); |
|
812 } |
|
813 _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end()); |
|
814 |
|
815 FilterMap filter_map(_res_cap); |
|
816 ResDigraph rgr(_sgr, filter_map); |
|
817 std::vector<int> cycle; |
|
818 std::vector<StaticDigraph::Arc> pred(_res_arc_num); |
|
819 PredMap pred_map(pred); |
|
820 BF bf(rgr, _cost_map); |
|
821 bf.distMap(_pi_map).predMap(pred_map); |
445 |
822 |
446 int length_bound = BF_FIRST_LIMIT; |
823 int length_bound = BF_FIRST_LIMIT; |
447 bool optimal = false; |
824 bool optimal = false; |
448 while (!optimal) { |
825 while (!optimal) { |
449 BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost); |
|
450 bf.predMap(pred); |
|
451 bf.init(0); |
826 bf.init(0); |
452 int iter_num = 0; |
827 int iter_num = 0; |
453 bool cycle_found = false; |
828 bool cycle_found = false; |
454 while (!cycle_found) { |
829 while (!cycle_found) { |
455 int curr_iter_num = iter_num + length_bound <= node_num ? |
830 // Perform some iterations of the Bellman-Ford algorithm |
456 length_bound : node_num - iter_num; |
831 int curr_iter_num = iter_num + length_bound <= _node_num ? |
|
832 length_bound : _node_num - iter_num; |
457 iter_num += curr_iter_num; |
833 iter_num += curr_iter_num; |
458 int real_iter_num = curr_iter_num; |
834 int real_iter_num = curr_iter_num; |
459 for (int i = 0; i < curr_iter_num; ++i) { |
835 for (int i = 0; i < curr_iter_num; ++i) { |
460 if (bf.processNextWeakRound()) { |
836 if (bf.processNextWeakRound()) { |
461 real_iter_num = i; |
837 real_iter_num = i; |
463 } |
839 } |
464 } |
840 } |
465 if (real_iter_num < curr_iter_num) { |
841 if (real_iter_num < curr_iter_num) { |
466 // Optimal flow is found |
842 // Optimal flow is found |
467 optimal = true; |
843 optimal = true; |
468 // Setting node potentials |
|
469 for (NodeIt n(_graph); n != INVALID; ++n) |
|
470 (*_potential)[n] = bf.dist(n); |
|
471 break; |
844 break; |
472 } else { |
845 } else { |
473 // Searching for node disjoint negative cycles |
846 // Search for node disjoint negative cycles |
474 for (ResNodeIt n(*_res_graph); n != INVALID; ++n) |
847 std::vector<int> state(_res_node_num, 0); |
475 visited[n] = 0; |
|
476 int id = 0; |
848 int id = 0; |
477 for (ResNodeIt n(*_res_graph); n != INVALID; ++n) { |
849 for (int u = 0; u != _res_node_num; ++u) { |
478 if (visited[n] > 0) continue; |
850 if (state[u] != 0) continue; |
479 visited[n] = ++id; |
851 ++id; |
480 ResNode u = pred[n] == INVALID ? |
852 int v = u; |
481 INVALID : _res_graph->source(pred[n]); |
853 for (; v != -1 && state[v] == 0; v = pred[v] == INVALID ? |
482 while (u != INVALID && visited[u] == 0) { |
854 -1 : rgr.id(rgr.source(pred[v]))) { |
483 visited[u] = id; |
855 state[v] = id; |
484 u = pred[u] == INVALID ? |
|
485 INVALID : _res_graph->source(pred[u]); |
|
486 } |
856 } |
487 if (u != INVALID && visited[u] == id) { |
857 if (v != -1 && state[v] == id) { |
488 // Finding the negative cycle |
858 // A negative cycle is found |
489 cycle_found = true; |
859 cycle_found = true; |
490 cycle.clear(); |
860 cycle.clear(); |
491 ResArc e = pred[u]; |
861 StaticDigraph::Arc a = pred[v]; |
492 cycle.push_back(e); |
862 Value d, delta = _res_cap[rgr.id(a)]; |
493 Capacity d = _res_graph->residualCapacity(e); |
863 cycle.push_back(rgr.id(a)); |
494 while (_res_graph->source(e) != u) { |
864 while (rgr.id(rgr.source(a)) != v) { |
495 cycle.push_back(e = pred[_res_graph->source(e)]); |
865 a = pred_map[rgr.source(a)]; |
496 if (_res_graph->residualCapacity(e) < d) |
866 d = _res_cap[rgr.id(a)]; |
497 d = _res_graph->residualCapacity(e); |
867 if (d < delta) delta = d; |
|
868 cycle.push_back(rgr.id(a)); |
498 } |
869 } |
499 |
870 |
500 // Augmenting along the cycle |
871 // Augment along the cycle |
501 for (int i = 0; i < int(cycle.size()); ++i) |
872 for (int i = 0; i < int(cycle.size()); ++i) { |
502 _res_graph->augment(cycle[i], d); |
873 int j = cycle[i]; |
|
874 _res_cap[j] -= delta; |
|
875 _res_cap[_reverse[j]] += delta; |
|
876 } |
503 } |
877 } |
504 } |
878 } |
505 } |
879 } |
506 |
880 |
507 if (!cycle_found) |
881 // Increase iteration limit if no cycle is found |
508 length_bound = length_bound * BF_LIMIT_FACTOR / 100; |
882 if (!cycle_found) { |
509 } |
883 length_bound = static_cast<int>(length_bound * BF_LIMIT_FACTOR); |
510 } |
884 } |
511 } |
885 } |
512 |
886 } |
513 /// \brief Execute the algorithm using \ref Howard. |
887 } |
514 /// |
888 |
515 /// Execute the algorithm using \ref Howard for negative |
889 // Execute the "Minimum Mean Cycle Canceling" method |
516 /// cycle detection. |
890 void startMinMeanCycleCanceling() { |
517 void startMinMean() { |
891 typedef SimplePath<StaticDigraph> SPath; |
518 typedef Path<ResDigraph> ResPath; |
892 typedef typename SPath::ArcIt SPathArcIt; |
519 Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost); |
893 typedef typename Howard<StaticDigraph, CostArcMap> |
520 ResPath cycle; |
894 ::template SetPath<SPath>::Create MMC; |
521 |
895 |
|
896 SPath cycle; |
|
897 MMC mmc(_sgr, _cost_map); |
522 mmc.cycle(cycle); |
898 mmc.cycle(cycle); |
523 if (mmc.findMinMean()) { |
899 buildResidualNetwork(); |
524 while (mmc.cycleLength() < 0) { |
900 while (mmc.findMinMean() && mmc.cycleLength() < 0) { |
525 // Finding the cycle |
901 // Find the cycle |
526 mmc.findCycle(); |
902 mmc.findCycle(); |
527 |
903 |
528 // Finding the largest flow amount that can be augmented |
904 // Compute delta value |
529 // along the cycle |
905 Value delta = INF; |
530 Capacity delta = 0; |
906 for (SPathArcIt a(cycle); a != INVALID; ++a) { |
531 for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) { |
907 Value d = _res_cap[_id_vec[_sgr.id(a)]]; |
532 if (delta == 0 || _res_graph->residualCapacity(e) < delta) |
908 if (d < delta) delta = d; |
533 delta = _res_graph->residualCapacity(e); |
909 } |
534 } |
910 |
535 |
911 // Augment along the cycle |
536 // Augmenting along the cycle |
912 for (SPathArcIt a(cycle); a != INVALID; ++a) { |
537 for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) |
913 int j = _id_vec[_sgr.id(a)]; |
538 _res_graph->augment(e, delta); |
914 _res_cap[j] -= delta; |
539 |
915 _res_cap[_reverse[j]] += delta; |
540 // Finding the minimum cycle mean for the modified residual |
916 } |
541 // digraph |
917 |
542 if (!mmc.findMinMean()) break; |
918 // Rebuild the residual network |
543 } |
919 buildResidualNetwork(); |
544 } |
920 } |
545 |
921 } |
546 // Computing node potentials |
922 |
547 BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost); |
923 // Execute the "Cancel And Tighten" method |
548 bf.init(0); bf.start(); |
924 void startCancelAndTighten() { |
549 for (NodeIt n(_graph); n != INVALID; ++n) |
925 // Constants for the min mean cycle computations |
550 (*_potential)[n] = bf.dist(n); |
926 const double LIMIT_FACTOR = 1.0; |
|
927 const int MIN_LIMIT = 5; |
|
928 |
|
929 // Contruct auxiliary data vectors |
|
930 DoubleVector pi(_res_node_num, 0.0); |
|
931 IntVector level(_res_node_num); |
|
932 CharVector reached(_res_node_num); |
|
933 CharVector processed(_res_node_num); |
|
934 IntVector pred_node(_res_node_num); |
|
935 IntVector pred_arc(_res_node_num); |
|
936 std::vector<int> stack(_res_node_num); |
|
937 std::vector<int> proc_vector(_res_node_num); |
|
938 |
|
939 // Initialize epsilon |
|
940 double epsilon = 0; |
|
941 for (int a = 0; a != _res_arc_num; ++a) { |
|
942 if (_res_cap[a] > 0 && -_cost[a] > epsilon) |
|
943 epsilon = -_cost[a]; |
|
944 } |
|
945 |
|
946 // Start phases |
|
947 Tolerance<double> tol; |
|
948 tol.epsilon(1e-6); |
|
949 int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num))); |
|
950 if (limit < MIN_LIMIT) limit = MIN_LIMIT; |
|
951 int iter = limit; |
|
952 while (epsilon * _res_node_num >= 1) { |
|
953 // Find and cancel cycles in the admissible network using DFS |
|
954 for (int u = 0; u != _res_node_num; ++u) { |
|
955 reached[u] = false; |
|
956 processed[u] = false; |
|
957 } |
|
958 int stack_head = -1; |
|
959 int proc_head = -1; |
|
960 for (int start = 0; start != _res_node_num; ++start) { |
|
961 if (reached[start]) continue; |
|
962 |
|
963 // New start node |
|
964 reached[start] = true; |
|
965 pred_arc[start] = -1; |
|
966 pred_node[start] = -1; |
|
967 |
|
968 // Find the first admissible outgoing arc |
|
969 double p = pi[start]; |
|
970 int a = _first_out[start]; |
|
971 int last_out = _first_out[start+1]; |
|
972 for (; a != last_out && (_res_cap[a] == 0 || |
|
973 !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ; |
|
974 if (a == last_out) { |
|
975 processed[start] = true; |
|
976 proc_vector[++proc_head] = start; |
|
977 continue; |
|
978 } |
|
979 stack[++stack_head] = a; |
|
980 |
|
981 while (stack_head >= 0) { |
|
982 int sa = stack[stack_head]; |
|
983 int u = _source[sa]; |
|
984 int v = _target[sa]; |
|
985 |
|
986 if (!reached[v]) { |
|
987 // A new node is reached |
|
988 reached[v] = true; |
|
989 pred_node[v] = u; |
|
990 pred_arc[v] = sa; |
|
991 p = pi[v]; |
|
992 a = _first_out[v]; |
|
993 last_out = _first_out[v+1]; |
|
994 for (; a != last_out && (_res_cap[a] == 0 || |
|
995 !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ; |
|
996 stack[++stack_head] = a == last_out ? -1 : a; |
|
997 } else { |
|
998 if (!processed[v]) { |
|
999 // A cycle is found |
|
1000 int n, w = u; |
|
1001 Value d, delta = _res_cap[sa]; |
|
1002 for (n = u; n != v; n = pred_node[n]) { |
|
1003 d = _res_cap[pred_arc[n]]; |
|
1004 if (d <= delta) { |
|
1005 delta = d; |
|
1006 w = pred_node[n]; |
|
1007 } |
|
1008 } |
|
1009 |
|
1010 // Augment along the cycle |
|
1011 _res_cap[sa] -= delta; |
|
1012 _res_cap[_reverse[sa]] += delta; |
|
1013 for (n = u; n != v; n = pred_node[n]) { |
|
1014 int pa = pred_arc[n]; |
|
1015 _res_cap[pa] -= delta; |
|
1016 _res_cap[_reverse[pa]] += delta; |
|
1017 } |
|
1018 for (n = u; stack_head > 0 && n != w; n = pred_node[n]) { |
|
1019 --stack_head; |
|
1020 reached[n] = false; |
|
1021 } |
|
1022 u = w; |
|
1023 } |
|
1024 v = u; |
|
1025 |
|
1026 // Find the next admissible outgoing arc |
|
1027 p = pi[v]; |
|
1028 a = stack[stack_head] + 1; |
|
1029 last_out = _first_out[v+1]; |
|
1030 for (; a != last_out && (_res_cap[a] == 0 || |
|
1031 !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ; |
|
1032 stack[stack_head] = a == last_out ? -1 : a; |
|
1033 } |
|
1034 |
|
1035 while (stack_head >= 0 && stack[stack_head] == -1) { |
|
1036 processed[v] = true; |
|
1037 proc_vector[++proc_head] = v; |
|
1038 if (--stack_head >= 0) { |
|
1039 // Find the next admissible outgoing arc |
|
1040 v = _source[stack[stack_head]]; |
|
1041 p = pi[v]; |
|
1042 a = stack[stack_head] + 1; |
|
1043 last_out = _first_out[v+1]; |
|
1044 for (; a != last_out && (_res_cap[a] == 0 || |
|
1045 !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ; |
|
1046 stack[stack_head] = a == last_out ? -1 : a; |
|
1047 } |
|
1048 } |
|
1049 } |
|
1050 } |
|
1051 |
|
1052 // Tighten potentials and epsilon |
|
1053 if (--iter > 0) { |
|
1054 for (int u = 0; u != _res_node_num; ++u) { |
|
1055 level[u] = 0; |
|
1056 } |
|
1057 for (int i = proc_head; i > 0; --i) { |
|
1058 int u = proc_vector[i]; |
|
1059 double p = pi[u]; |
|
1060 int l = level[u] + 1; |
|
1061 int last_out = _first_out[u+1]; |
|
1062 for (int a = _first_out[u]; a != last_out; ++a) { |
|
1063 int v = _target[a]; |
|
1064 if (_res_cap[a] > 0 && tol.negative(_cost[a] + p - pi[v]) && |
|
1065 l > level[v]) level[v] = l; |
|
1066 } |
|
1067 } |
|
1068 |
|
1069 // Modify potentials |
|
1070 double q = std::numeric_limits<double>::max(); |
|
1071 for (int u = 0; u != _res_node_num; ++u) { |
|
1072 int lu = level[u]; |
|
1073 double p, pu = pi[u]; |
|
1074 int last_out = _first_out[u+1]; |
|
1075 for (int a = _first_out[u]; a != last_out; ++a) { |
|
1076 if (_res_cap[a] == 0) continue; |
|
1077 int v = _target[a]; |
|
1078 int ld = lu - level[v]; |
|
1079 if (ld > 0) { |
|
1080 p = (_cost[a] + pu - pi[v] + epsilon) / (ld + 1); |
|
1081 if (p < q) q = p; |
|
1082 } |
|
1083 } |
|
1084 } |
|
1085 for (int u = 0; u != _res_node_num; ++u) { |
|
1086 pi[u] -= q * level[u]; |
|
1087 } |
|
1088 |
|
1089 // Modify epsilon |
|
1090 epsilon = 0; |
|
1091 for (int u = 0; u != _res_node_num; ++u) { |
|
1092 double curr, pu = pi[u]; |
|
1093 int last_out = _first_out[u+1]; |
|
1094 for (int a = _first_out[u]; a != last_out; ++a) { |
|
1095 if (_res_cap[a] == 0) continue; |
|
1096 curr = _cost[a] + pu - pi[_target[a]]; |
|
1097 if (-curr > epsilon) epsilon = -curr; |
|
1098 } |
|
1099 } |
|
1100 } else { |
|
1101 typedef Howard<StaticDigraph, CostArcMap> MMC; |
|
1102 typedef typename BellmanFord<StaticDigraph, CostArcMap> |
|
1103 ::template SetDistMap<CostNodeMap>::Create BF; |
|
1104 |
|
1105 // Set epsilon to the minimum cycle mean |
|
1106 buildResidualNetwork(); |
|
1107 MMC mmc(_sgr, _cost_map); |
|
1108 mmc.findMinMean(); |
|
1109 epsilon = -mmc.cycleMean(); |
|
1110 Cost cycle_cost = mmc.cycleLength(); |
|
1111 int cycle_size = mmc.cycleArcNum(); |
|
1112 |
|
1113 // Compute feasible potentials for the current epsilon |
|
1114 for (int i = 0; i != int(_cost_vec.size()); ++i) { |
|
1115 _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost; |
|
1116 } |
|
1117 BF bf(_sgr, _cost_map); |
|
1118 bf.distMap(_pi_map); |
|
1119 bf.init(0); |
|
1120 bf.start(); |
|
1121 for (int u = 0; u != _res_node_num; ++u) { |
|
1122 pi[u] = static_cast<double>(_pi[u]) / cycle_size; |
|
1123 } |
|
1124 |
|
1125 iter = limit; |
|
1126 } |
|
1127 } |
551 } |
1128 } |
552 |
1129 |
553 }; //class CycleCanceling |
1130 }; //class CycleCanceling |
554 |
1131 |
555 ///@} |
1132 ///@} |