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