kpeter@808
|
1 |
/* -*- C++ -*-
|
kpeter@808
|
2 |
*
|
kpeter@808
|
3 |
* This file is a part of LEMON, a generic C++ optimization library
|
kpeter@808
|
4 |
*
|
kpeter@808
|
5 |
* Copyright (C) 2003-2008
|
kpeter@808
|
6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
|
kpeter@808
|
7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES).
|
kpeter@808
|
8 |
*
|
kpeter@808
|
9 |
* Permission to use, modify and distribute this software is granted
|
kpeter@808
|
10 |
* provided that this copyright notice appears in all copies. For
|
kpeter@808
|
11 |
* precise terms see the accompanying LICENSE file.
|
kpeter@808
|
12 |
*
|
kpeter@808
|
13 |
* This software is provided "AS IS" with no warranty of any kind,
|
kpeter@808
|
14 |
* express or implied, and with no claim as to its suitability for any
|
kpeter@808
|
15 |
* purpose.
|
kpeter@808
|
16 |
*
|
kpeter@808
|
17 |
*/
|
kpeter@808
|
18 |
|
kpeter@808
|
19 |
#ifndef LEMON_COST_SCALING_H
|
kpeter@808
|
20 |
#define LEMON_COST_SCALING_H
|
kpeter@808
|
21 |
|
kpeter@808
|
22 |
/// \ingroup min_cost_flow_algs
|
kpeter@808
|
23 |
/// \file
|
kpeter@808
|
24 |
/// \brief Cost scaling algorithm for finding a minimum cost flow.
|
kpeter@808
|
25 |
|
kpeter@808
|
26 |
#include <vector>
|
kpeter@808
|
27 |
#include <deque>
|
kpeter@808
|
28 |
#include <limits>
|
kpeter@808
|
29 |
|
kpeter@808
|
30 |
#include <lemon/core.h>
|
kpeter@808
|
31 |
#include <lemon/maps.h>
|
kpeter@808
|
32 |
#include <lemon/math.h>
|
kpeter@808
|
33 |
#include <lemon/adaptors.h>
|
kpeter@808
|
34 |
#include <lemon/circulation.h>
|
kpeter@808
|
35 |
#include <lemon/bellman_ford.h>
|
kpeter@808
|
36 |
|
kpeter@808
|
37 |
namespace lemon {
|
kpeter@808
|
38 |
|
kpeter@808
|
39 |
/// \addtogroup min_cost_flow_algs
|
kpeter@808
|
40 |
/// @{
|
kpeter@808
|
41 |
|
kpeter@808
|
42 |
/// \brief Implementation of the cost scaling algorithm for finding a
|
kpeter@808
|
43 |
/// minimum cost flow.
|
kpeter@808
|
44 |
///
|
kpeter@808
|
45 |
/// \ref CostScaling implements the cost scaling algorithm performing
|
kpeter@808
|
46 |
/// augment/push and relabel operations for finding a minimum cost
|
kpeter@808
|
47 |
/// flow.
|
kpeter@808
|
48 |
///
|
kpeter@808
|
49 |
/// \tparam Digraph The digraph type the algorithm runs on.
|
kpeter@808
|
50 |
/// \tparam LowerMap The type of the lower bound map.
|
kpeter@808
|
51 |
/// \tparam CapacityMap The type of the capacity (upper bound) map.
|
kpeter@808
|
52 |
/// \tparam CostMap The type of the cost (length) map.
|
kpeter@808
|
53 |
/// \tparam SupplyMap The type of the supply map.
|
kpeter@808
|
54 |
///
|
kpeter@808
|
55 |
/// \warning
|
kpeter@808
|
56 |
/// - Arc capacities and costs should be \e non-negative \e integers.
|
kpeter@808
|
57 |
/// - Supply values should be \e signed \e integers.
|
kpeter@808
|
58 |
/// - The value types of the maps should be convertible to each other.
|
kpeter@808
|
59 |
/// - \c CostMap::Value must be signed type.
|
kpeter@808
|
60 |
///
|
kpeter@808
|
61 |
/// \note Arc costs are multiplied with the number of nodes during
|
kpeter@808
|
62 |
/// the algorithm so overflow problems may arise more easily than with
|
kpeter@808
|
63 |
/// other minimum cost flow algorithms.
|
kpeter@808
|
64 |
/// If it is available, <tt>long long int</tt> type is used instead of
|
kpeter@808
|
65 |
/// <tt>long int</tt> in the inside computations.
|
kpeter@808
|
66 |
///
|
kpeter@808
|
67 |
/// \author Peter Kovacs
|
kpeter@808
|
68 |
template < typename Digraph,
|
kpeter@808
|
69 |
typename LowerMap = typename Digraph::template ArcMap<int>,
|
kpeter@808
|
70 |
typename CapacityMap = typename Digraph::template ArcMap<int>,
|
kpeter@808
|
71 |
typename CostMap = typename Digraph::template ArcMap<int>,
|
kpeter@808
|
72 |
typename SupplyMap = typename Digraph::template NodeMap<int> >
|
kpeter@808
|
73 |
class CostScaling
|
kpeter@808
|
74 |
{
|
kpeter@808
|
75 |
TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
|
kpeter@808
|
76 |
|
kpeter@808
|
77 |
typedef typename CapacityMap::Value Capacity;
|
kpeter@808
|
78 |
typedef typename CostMap::Value Cost;
|
kpeter@808
|
79 |
typedef typename SupplyMap::Value Supply;
|
kpeter@808
|
80 |
typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
|
kpeter@808
|
81 |
typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
|
kpeter@808
|
82 |
|
kpeter@808
|
83 |
typedef ResidualDigraph< const Digraph,
|
kpeter@808
|
84 |
CapacityArcMap, CapacityArcMap > ResDigraph;
|
kpeter@808
|
85 |
typedef typename ResDigraph::Arc ResArc;
|
kpeter@808
|
86 |
|
kpeter@808
|
87 |
#if defined __GNUC__ && !defined __STRICT_ANSI__
|
kpeter@808
|
88 |
typedef long long int LCost;
|
kpeter@808
|
89 |
#else
|
kpeter@808
|
90 |
typedef long int LCost;
|
kpeter@808
|
91 |
#endif
|
kpeter@808
|
92 |
typedef typename Digraph::template ArcMap<LCost> LargeCostMap;
|
kpeter@808
|
93 |
|
kpeter@808
|
94 |
public:
|
kpeter@808
|
95 |
|
kpeter@808
|
96 |
/// The type of the flow map.
|
kpeter@808
|
97 |
typedef typename Digraph::template ArcMap<Capacity> FlowMap;
|
kpeter@808
|
98 |
/// The type of the potential map.
|
kpeter@808
|
99 |
typedef typename Digraph::template NodeMap<LCost> PotentialMap;
|
kpeter@808
|
100 |
|
kpeter@808
|
101 |
private:
|
kpeter@808
|
102 |
|
kpeter@808
|
103 |
/// \brief Map adaptor class for handling residual arc costs.
|
kpeter@808
|
104 |
///
|
kpeter@808
|
105 |
/// Map adaptor class for handling residual arc costs.
|
kpeter@808
|
106 |
template <typename Map>
|
kpeter@808
|
107 |
class ResidualCostMap : public MapBase<ResArc, typename Map::Value>
|
kpeter@808
|
108 |
{
|
kpeter@808
|
109 |
private:
|
kpeter@808
|
110 |
|
kpeter@808
|
111 |
const Map &_cost_map;
|
kpeter@808
|
112 |
|
kpeter@808
|
113 |
public:
|
kpeter@808
|
114 |
|
kpeter@808
|
115 |
///\e
|
kpeter@808
|
116 |
ResidualCostMap(const Map &cost_map) :
|
kpeter@808
|
117 |
_cost_map(cost_map) {}
|
kpeter@808
|
118 |
|
kpeter@808
|
119 |
///\e
|
kpeter@808
|
120 |
inline typename Map::Value operator[](const ResArc &e) const {
|
kpeter@808
|
121 |
return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
|
kpeter@808
|
122 |
}
|
kpeter@808
|
123 |
|
kpeter@808
|
124 |
}; //class ResidualCostMap
|
kpeter@808
|
125 |
|
kpeter@808
|
126 |
/// \brief Map adaptor class for handling reduced arc costs.
|
kpeter@808
|
127 |
///
|
kpeter@808
|
128 |
/// Map adaptor class for handling reduced arc costs.
|
kpeter@808
|
129 |
class ReducedCostMap : public MapBase<Arc, LCost>
|
kpeter@808
|
130 |
{
|
kpeter@808
|
131 |
private:
|
kpeter@808
|
132 |
|
kpeter@808
|
133 |
const Digraph &_gr;
|
kpeter@808
|
134 |
const LargeCostMap &_cost_map;
|
kpeter@808
|
135 |
const PotentialMap &_pot_map;
|
kpeter@808
|
136 |
|
kpeter@808
|
137 |
public:
|
kpeter@808
|
138 |
|
kpeter@808
|
139 |
///\e
|
kpeter@808
|
140 |
ReducedCostMap( const Digraph &gr,
|
kpeter@808
|
141 |
const LargeCostMap &cost_map,
|
kpeter@808
|
142 |
const PotentialMap &pot_map ) :
|
kpeter@808
|
143 |
_gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
|
kpeter@808
|
144 |
|
kpeter@808
|
145 |
///\e
|
kpeter@808
|
146 |
inline LCost operator[](const Arc &e) const {
|
kpeter@808
|
147 |
return _cost_map[e] + _pot_map[_gr.source(e)]
|
kpeter@808
|
148 |
- _pot_map[_gr.target(e)];
|
kpeter@808
|
149 |
}
|
kpeter@808
|
150 |
|
kpeter@808
|
151 |
}; //class ReducedCostMap
|
kpeter@808
|
152 |
|
kpeter@808
|
153 |
private:
|
kpeter@808
|
154 |
|
kpeter@808
|
155 |
// The digraph the algorithm runs on
|
kpeter@808
|
156 |
const Digraph &_graph;
|
kpeter@808
|
157 |
// The original lower bound map
|
kpeter@808
|
158 |
const LowerMap *_lower;
|
kpeter@808
|
159 |
// The modified capacity map
|
kpeter@808
|
160 |
CapacityArcMap _capacity;
|
kpeter@808
|
161 |
// The original cost map
|
kpeter@808
|
162 |
const CostMap &_orig_cost;
|
kpeter@808
|
163 |
// The scaled cost map
|
kpeter@808
|
164 |
LargeCostMap _cost;
|
kpeter@808
|
165 |
// The modified supply map
|
kpeter@808
|
166 |
SupplyNodeMap _supply;
|
kpeter@808
|
167 |
bool _valid_supply;
|
kpeter@808
|
168 |
|
kpeter@808
|
169 |
// Arc map of the current flow
|
kpeter@808
|
170 |
FlowMap *_flow;
|
kpeter@808
|
171 |
bool _local_flow;
|
kpeter@808
|
172 |
// Node map of the current potentials
|
kpeter@808
|
173 |
PotentialMap *_potential;
|
kpeter@808
|
174 |
bool _local_potential;
|
kpeter@808
|
175 |
|
kpeter@808
|
176 |
// The residual cost map
|
kpeter@808
|
177 |
ResidualCostMap<LargeCostMap> _res_cost;
|
kpeter@808
|
178 |
// The residual digraph
|
kpeter@808
|
179 |
ResDigraph *_res_graph;
|
kpeter@808
|
180 |
// The reduced cost map
|
kpeter@808
|
181 |
ReducedCostMap *_red_cost;
|
kpeter@808
|
182 |
// The excess map
|
kpeter@808
|
183 |
SupplyNodeMap _excess;
|
kpeter@808
|
184 |
// The epsilon parameter used for cost scaling
|
kpeter@808
|
185 |
LCost _epsilon;
|
kpeter@808
|
186 |
// The scaling factor
|
kpeter@808
|
187 |
int _alpha;
|
kpeter@808
|
188 |
|
kpeter@808
|
189 |
public:
|
kpeter@808
|
190 |
|
kpeter@808
|
191 |
/// \brief General constructor (with lower bounds).
|
kpeter@808
|
192 |
///
|
kpeter@808
|
193 |
/// General constructor (with lower bounds).
|
kpeter@808
|
194 |
///
|
kpeter@808
|
195 |
/// \param digraph The digraph the algorithm runs on.
|
kpeter@808
|
196 |
/// \param lower The lower bounds of the arcs.
|
kpeter@808
|
197 |
/// \param capacity The capacities (upper bounds) of the arcs.
|
kpeter@808
|
198 |
/// \param cost The cost (length) values of the arcs.
|
kpeter@808
|
199 |
/// \param supply The supply values of the nodes (signed).
|
kpeter@808
|
200 |
CostScaling( const Digraph &digraph,
|
kpeter@808
|
201 |
const LowerMap &lower,
|
kpeter@808
|
202 |
const CapacityMap &capacity,
|
kpeter@808
|
203 |
const CostMap &cost,
|
kpeter@808
|
204 |
const SupplyMap &supply ) :
|
kpeter@808
|
205 |
_graph(digraph), _lower(&lower), _capacity(digraph), _orig_cost(cost),
|
kpeter@808
|
206 |
_cost(digraph), _supply(digraph), _flow(NULL), _local_flow(false),
|
kpeter@808
|
207 |
_potential(NULL), _local_potential(false), _res_cost(_cost),
|
kpeter@808
|
208 |
_res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
|
kpeter@808
|
209 |
{
|
kpeter@808
|
210 |
// Check the sum of supply values
|
kpeter@808
|
211 |
Supply sum = 0;
|
kpeter@808
|
212 |
for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
|
kpeter@808
|
213 |
_valid_supply = sum == 0;
|
kpeter@808
|
214 |
|
kpeter@808
|
215 |
for (ArcIt e(_graph); e != INVALID; ++e) _capacity[e] = capacity[e];
|
kpeter@808
|
216 |
for (NodeIt n(_graph); n != INVALID; ++n) _supply[n] = supply[n];
|
kpeter@808
|
217 |
|
kpeter@808
|
218 |
// Remove non-zero lower bounds
|
kpeter@808
|
219 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
kpeter@808
|
220 |
if (lower[e] != 0) {
|
kpeter@808
|
221 |
_capacity[e] -= lower[e];
|
kpeter@808
|
222 |
_supply[_graph.source(e)] -= lower[e];
|
kpeter@808
|
223 |
_supply[_graph.target(e)] += lower[e];
|
kpeter@808
|
224 |
}
|
kpeter@808
|
225 |
}
|
kpeter@808
|
226 |
}
|
kpeter@808
|
227 |
/*
|
kpeter@808
|
228 |
/// \brief General constructor (without lower bounds).
|
kpeter@808
|
229 |
///
|
kpeter@808
|
230 |
/// General constructor (without lower bounds).
|
kpeter@808
|
231 |
///
|
kpeter@808
|
232 |
/// \param digraph The digraph the algorithm runs on.
|
kpeter@808
|
233 |
/// \param capacity The capacities (upper bounds) of the arcs.
|
kpeter@808
|
234 |
/// \param cost The cost (length) values of the arcs.
|
kpeter@808
|
235 |
/// \param supply The supply values of the nodes (signed).
|
kpeter@808
|
236 |
CostScaling( const Digraph &digraph,
|
kpeter@808
|
237 |
const CapacityMap &capacity,
|
kpeter@808
|
238 |
const CostMap &cost,
|
kpeter@808
|
239 |
const SupplyMap &supply ) :
|
kpeter@808
|
240 |
_graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
|
kpeter@808
|
241 |
_cost(digraph), _supply(supply), _flow(NULL), _local_flow(false),
|
kpeter@808
|
242 |
_potential(NULL), _local_potential(false), _res_cost(_cost),
|
kpeter@808
|
243 |
_res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
|
kpeter@808
|
244 |
{
|
kpeter@808
|
245 |
// Check the sum of supply values
|
kpeter@808
|
246 |
Supply sum = 0;
|
kpeter@808
|
247 |
for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
|
kpeter@808
|
248 |
_valid_supply = sum == 0;
|
kpeter@808
|
249 |
}
|
kpeter@808
|
250 |
|
kpeter@808
|
251 |
/// \brief Simple constructor (with lower bounds).
|
kpeter@808
|
252 |
///
|
kpeter@808
|
253 |
/// Simple constructor (with lower bounds).
|
kpeter@808
|
254 |
///
|
kpeter@808
|
255 |
/// \param digraph The digraph the algorithm runs on.
|
kpeter@808
|
256 |
/// \param lower The lower bounds of the arcs.
|
kpeter@808
|
257 |
/// \param capacity The capacities (upper bounds) of the arcs.
|
kpeter@808
|
258 |
/// \param cost The cost (length) values of the arcs.
|
kpeter@808
|
259 |
/// \param s The source node.
|
kpeter@808
|
260 |
/// \param t The target node.
|
kpeter@808
|
261 |
/// \param flow_value The required amount of flow from node \c s
|
kpeter@808
|
262 |
/// to node \c t (i.e. the supply of \c s and the demand of \c t).
|
kpeter@808
|
263 |
CostScaling( const Digraph &digraph,
|
kpeter@808
|
264 |
const LowerMap &lower,
|
kpeter@808
|
265 |
const CapacityMap &capacity,
|
kpeter@808
|
266 |
const CostMap &cost,
|
kpeter@808
|
267 |
Node s, Node t,
|
kpeter@808
|
268 |
Supply flow_value ) :
|
kpeter@808
|
269 |
_graph(digraph), _lower(&lower), _capacity(capacity), _orig_cost(cost),
|
kpeter@808
|
270 |
_cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
|
kpeter@808
|
271 |
_potential(NULL), _local_potential(false), _res_cost(_cost),
|
kpeter@808
|
272 |
_res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
|
kpeter@808
|
273 |
{
|
kpeter@808
|
274 |
// Remove non-zero lower bounds
|
kpeter@808
|
275 |
_supply[s] = flow_value;
|
kpeter@808
|
276 |
_supply[t] = -flow_value;
|
kpeter@808
|
277 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
kpeter@808
|
278 |
if (lower[e] != 0) {
|
kpeter@808
|
279 |
_capacity[e] -= lower[e];
|
kpeter@808
|
280 |
_supply[_graph.source(e)] -= lower[e];
|
kpeter@808
|
281 |
_supply[_graph.target(e)] += lower[e];
|
kpeter@808
|
282 |
}
|
kpeter@808
|
283 |
}
|
kpeter@808
|
284 |
_valid_supply = true;
|
kpeter@808
|
285 |
}
|
kpeter@808
|
286 |
|
kpeter@808
|
287 |
/// \brief Simple constructor (without lower bounds).
|
kpeter@808
|
288 |
///
|
kpeter@808
|
289 |
/// Simple constructor (without lower bounds).
|
kpeter@808
|
290 |
///
|
kpeter@808
|
291 |
/// \param digraph The digraph the algorithm runs on.
|
kpeter@808
|
292 |
/// \param capacity The capacities (upper bounds) of the arcs.
|
kpeter@808
|
293 |
/// \param cost The cost (length) values of the arcs.
|
kpeter@808
|
294 |
/// \param s The source node.
|
kpeter@808
|
295 |
/// \param t The target node.
|
kpeter@808
|
296 |
/// \param flow_value The required amount of flow from node \c s
|
kpeter@808
|
297 |
/// to node \c t (i.e. the supply of \c s and the demand of \c t).
|
kpeter@808
|
298 |
CostScaling( const Digraph &digraph,
|
kpeter@808
|
299 |
const CapacityMap &capacity,
|
kpeter@808
|
300 |
const CostMap &cost,
|
kpeter@808
|
301 |
Node s, Node t,
|
kpeter@808
|
302 |
Supply flow_value ) :
|
kpeter@808
|
303 |
_graph(digraph), _lower(NULL), _capacity(capacity), _orig_cost(cost),
|
kpeter@808
|
304 |
_cost(digraph), _supply(digraph, 0), _flow(NULL), _local_flow(false),
|
kpeter@808
|
305 |
_potential(NULL), _local_potential(false), _res_cost(_cost),
|
kpeter@808
|
306 |
_res_graph(NULL), _red_cost(NULL), _excess(digraph, 0)
|
kpeter@808
|
307 |
{
|
kpeter@808
|
308 |
_supply[s] = flow_value;
|
kpeter@808
|
309 |
_supply[t] = -flow_value;
|
kpeter@808
|
310 |
_valid_supply = true;
|
kpeter@808
|
311 |
}
|
kpeter@808
|
312 |
*/
|
kpeter@808
|
313 |
/// Destructor.
|
kpeter@808
|
314 |
~CostScaling() {
|
kpeter@808
|
315 |
if (_local_flow) delete _flow;
|
kpeter@808
|
316 |
if (_local_potential) delete _potential;
|
kpeter@808
|
317 |
delete _res_graph;
|
kpeter@808
|
318 |
delete _red_cost;
|
kpeter@808
|
319 |
}
|
kpeter@808
|
320 |
|
kpeter@808
|
321 |
/// \brief Set the flow map.
|
kpeter@808
|
322 |
///
|
kpeter@808
|
323 |
/// Set the flow map.
|
kpeter@808
|
324 |
///
|
kpeter@808
|
325 |
/// \return \c (*this)
|
kpeter@808
|
326 |
CostScaling& flowMap(FlowMap &map) {
|
kpeter@808
|
327 |
if (_local_flow) {
|
kpeter@808
|
328 |
delete _flow;
|
kpeter@808
|
329 |
_local_flow = false;
|
kpeter@808
|
330 |
}
|
kpeter@808
|
331 |
_flow = ↦
|
kpeter@808
|
332 |
return *this;
|
kpeter@808
|
333 |
}
|
kpeter@808
|
334 |
|
kpeter@808
|
335 |
/// \brief Set the potential map.
|
kpeter@808
|
336 |
///
|
kpeter@808
|
337 |
/// Set the potential map.
|
kpeter@808
|
338 |
///
|
kpeter@808
|
339 |
/// \return \c (*this)
|
kpeter@808
|
340 |
CostScaling& potentialMap(PotentialMap &map) {
|
kpeter@808
|
341 |
if (_local_potential) {
|
kpeter@808
|
342 |
delete _potential;
|
kpeter@808
|
343 |
_local_potential = false;
|
kpeter@808
|
344 |
}
|
kpeter@808
|
345 |
_potential = ↦
|
kpeter@808
|
346 |
return *this;
|
kpeter@808
|
347 |
}
|
kpeter@808
|
348 |
|
kpeter@808
|
349 |
/// \name Execution control
|
kpeter@808
|
350 |
|
kpeter@808
|
351 |
/// @{
|
kpeter@808
|
352 |
|
kpeter@808
|
353 |
/// \brief Run the algorithm.
|
kpeter@808
|
354 |
///
|
kpeter@808
|
355 |
/// Run the algorithm.
|
kpeter@808
|
356 |
///
|
kpeter@808
|
357 |
/// \param partial_augment By default the algorithm performs
|
kpeter@808
|
358 |
/// partial augment and relabel operations in the cost scaling
|
kpeter@808
|
359 |
/// phases. Set this parameter to \c false for using local push and
|
kpeter@808
|
360 |
/// relabel operations instead.
|
kpeter@808
|
361 |
///
|
kpeter@808
|
362 |
/// \return \c true if a feasible flow can be found.
|
kpeter@808
|
363 |
bool run(bool partial_augment = true) {
|
kpeter@808
|
364 |
if (partial_augment) {
|
kpeter@808
|
365 |
return init() && startPartialAugment();
|
kpeter@808
|
366 |
} else {
|
kpeter@808
|
367 |
return init() && startPushRelabel();
|
kpeter@808
|
368 |
}
|
kpeter@808
|
369 |
}
|
kpeter@808
|
370 |
|
kpeter@808
|
371 |
/// @}
|
kpeter@808
|
372 |
|
kpeter@808
|
373 |
/// \name Query Functions
|
kpeter@808
|
374 |
/// The result of the algorithm can be obtained using these
|
kpeter@808
|
375 |
/// functions.\n
|
kpeter@808
|
376 |
/// \ref lemon::CostScaling::run() "run()" must be called before
|
kpeter@808
|
377 |
/// using them.
|
kpeter@808
|
378 |
|
kpeter@808
|
379 |
/// @{
|
kpeter@808
|
380 |
|
kpeter@808
|
381 |
/// \brief Return a const reference to the arc map storing the
|
kpeter@808
|
382 |
/// found flow.
|
kpeter@808
|
383 |
///
|
kpeter@808
|
384 |
/// Return a const reference to the arc map storing the found flow.
|
kpeter@808
|
385 |
///
|
kpeter@808
|
386 |
/// \pre \ref run() must be called before using this function.
|
kpeter@808
|
387 |
const FlowMap& flowMap() const {
|
kpeter@808
|
388 |
return *_flow;
|
kpeter@808
|
389 |
}
|
kpeter@808
|
390 |
|
kpeter@808
|
391 |
/// \brief Return a const reference to the node map storing the
|
kpeter@808
|
392 |
/// found potentials (the dual solution).
|
kpeter@808
|
393 |
///
|
kpeter@808
|
394 |
/// Return a const reference to the node map storing the found
|
kpeter@808
|
395 |
/// potentials (the dual solution).
|
kpeter@808
|
396 |
///
|
kpeter@808
|
397 |
/// \pre \ref run() must be called before using this function.
|
kpeter@808
|
398 |
const PotentialMap& potentialMap() const {
|
kpeter@808
|
399 |
return *_potential;
|
kpeter@808
|
400 |
}
|
kpeter@808
|
401 |
|
kpeter@808
|
402 |
/// \brief Return the flow on the given arc.
|
kpeter@808
|
403 |
///
|
kpeter@808
|
404 |
/// Return the flow on the given arc.
|
kpeter@808
|
405 |
///
|
kpeter@808
|
406 |
/// \pre \ref run() must be called before using this function.
|
kpeter@808
|
407 |
Capacity flow(const Arc& arc) const {
|
kpeter@808
|
408 |
return (*_flow)[arc];
|
kpeter@808
|
409 |
}
|
kpeter@808
|
410 |
|
kpeter@808
|
411 |
/// \brief Return the potential of the given node.
|
kpeter@808
|
412 |
///
|
kpeter@808
|
413 |
/// Return the potential of the given node.
|
kpeter@808
|
414 |
///
|
kpeter@808
|
415 |
/// \pre \ref run() must be called before using this function.
|
kpeter@808
|
416 |
Cost potential(const Node& node) const {
|
kpeter@808
|
417 |
return (*_potential)[node];
|
kpeter@808
|
418 |
}
|
kpeter@808
|
419 |
|
kpeter@808
|
420 |
/// \brief Return the total cost of the found flow.
|
kpeter@808
|
421 |
///
|
kpeter@808
|
422 |
/// Return the total cost of the found flow. The complexity of the
|
kpeter@808
|
423 |
/// function is \f$ O(e) \f$.
|
kpeter@808
|
424 |
///
|
kpeter@808
|
425 |
/// \pre \ref run() must be called before using this function.
|
kpeter@808
|
426 |
Cost totalCost() const {
|
kpeter@808
|
427 |
Cost c = 0;
|
kpeter@808
|
428 |
for (ArcIt e(_graph); e != INVALID; ++e)
|
kpeter@808
|
429 |
c += (*_flow)[e] * _orig_cost[e];
|
kpeter@808
|
430 |
return c;
|
kpeter@808
|
431 |
}
|
kpeter@808
|
432 |
|
kpeter@808
|
433 |
/// @}
|
kpeter@808
|
434 |
|
kpeter@808
|
435 |
private:
|
kpeter@808
|
436 |
|
kpeter@808
|
437 |
/// Initialize the algorithm.
|
kpeter@808
|
438 |
bool init() {
|
kpeter@808
|
439 |
if (!_valid_supply) return false;
|
kpeter@808
|
440 |
// The scaling factor
|
kpeter@808
|
441 |
_alpha = 8;
|
kpeter@808
|
442 |
|
kpeter@808
|
443 |
// Initialize flow and potential maps
|
kpeter@808
|
444 |
if (!_flow) {
|
kpeter@808
|
445 |
_flow = new FlowMap(_graph);
|
kpeter@808
|
446 |
_local_flow = true;
|
kpeter@808
|
447 |
}
|
kpeter@808
|
448 |
if (!_potential) {
|
kpeter@808
|
449 |
_potential = new PotentialMap(_graph);
|
kpeter@808
|
450 |
_local_potential = true;
|
kpeter@808
|
451 |
}
|
kpeter@808
|
452 |
|
kpeter@808
|
453 |
_red_cost = new ReducedCostMap(_graph, _cost, *_potential);
|
kpeter@808
|
454 |
_res_graph = new ResDigraph(_graph, _capacity, *_flow);
|
kpeter@808
|
455 |
|
kpeter@808
|
456 |
// Initialize the scaled cost map and the epsilon parameter
|
kpeter@808
|
457 |
Cost max_cost = 0;
|
kpeter@808
|
458 |
int node_num = countNodes(_graph);
|
kpeter@808
|
459 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
kpeter@808
|
460 |
_cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
|
kpeter@808
|
461 |
if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
|
kpeter@808
|
462 |
}
|
kpeter@808
|
463 |
_epsilon = max_cost * node_num;
|
kpeter@808
|
464 |
|
kpeter@808
|
465 |
// Find a feasible flow using Circulation
|
kpeter@808
|
466 |
Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
|
kpeter@808
|
467 |
SupplyMap >
|
kpeter@808
|
468 |
circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
|
kpeter@808
|
469 |
_supply );
|
kpeter@808
|
470 |
return circulation.flowMap(*_flow).run();
|
kpeter@808
|
471 |
}
|
kpeter@808
|
472 |
|
kpeter@808
|
473 |
/// Execute the algorithm performing partial augmentation and
|
kpeter@808
|
474 |
/// relabel operations.
|
kpeter@808
|
475 |
bool startPartialAugment() {
|
kpeter@808
|
476 |
// Paramters for heuristics
|
kpeter@808
|
477 |
// const int BF_HEURISTIC_EPSILON_BOUND = 1000;
|
kpeter@808
|
478 |
// const int BF_HEURISTIC_BOUND_FACTOR = 3;
|
kpeter@808
|
479 |
// Maximum augment path length
|
kpeter@808
|
480 |
const int MAX_PATH_LENGTH = 4;
|
kpeter@808
|
481 |
|
kpeter@808
|
482 |
// Variables
|
kpeter@808
|
483 |
typename Digraph::template NodeMap<Arc> pred_arc(_graph);
|
kpeter@808
|
484 |
typename Digraph::template NodeMap<bool> forward(_graph);
|
kpeter@808
|
485 |
typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
|
kpeter@808
|
486 |
typename Digraph::template NodeMap<InArcIt> next_in(_graph);
|
kpeter@808
|
487 |
typename Digraph::template NodeMap<bool> next_dir(_graph);
|
kpeter@808
|
488 |
std::deque<Node> active_nodes;
|
kpeter@808
|
489 |
std::vector<Node> path_nodes;
|
kpeter@808
|
490 |
|
kpeter@808
|
491 |
// int node_num = countNodes(_graph);
|
kpeter@808
|
492 |
for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
|
kpeter@808
|
493 |
1 : _epsilon / _alpha )
|
kpeter@808
|
494 |
{
|
kpeter@808
|
495 |
/*
|
kpeter@808
|
496 |
// "Early Termination" heuristic: use Bellman-Ford algorithm
|
kpeter@808
|
497 |
// to check if the current flow is optimal
|
kpeter@808
|
498 |
if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
|
kpeter@808
|
499 |
typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
|
kpeter@808
|
500 |
ShiftCostMap shift_cost(_res_cost, 1);
|
kpeter@808
|
501 |
BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
|
kpeter@808
|
502 |
bf.init(0);
|
kpeter@808
|
503 |
bool done = false;
|
kpeter@808
|
504 |
int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
|
kpeter@808
|
505 |
for (int i = 0; i < K && !done; ++i)
|
kpeter@808
|
506 |
done = bf.processNextWeakRound();
|
kpeter@808
|
507 |
if (done) break;
|
kpeter@808
|
508 |
}
|
kpeter@808
|
509 |
*/
|
kpeter@808
|
510 |
// Saturate arcs not satisfying the optimality condition
|
kpeter@808
|
511 |
Capacity delta;
|
kpeter@808
|
512 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
kpeter@808
|
513 |
if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
|
kpeter@808
|
514 |
delta = _capacity[e] - (*_flow)[e];
|
kpeter@808
|
515 |
_excess[_graph.source(e)] -= delta;
|
kpeter@808
|
516 |
_excess[_graph.target(e)] += delta;
|
kpeter@808
|
517 |
(*_flow)[e] = _capacity[e];
|
kpeter@808
|
518 |
}
|
kpeter@808
|
519 |
if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
|
kpeter@808
|
520 |
_excess[_graph.target(e)] -= (*_flow)[e];
|
kpeter@808
|
521 |
_excess[_graph.source(e)] += (*_flow)[e];
|
kpeter@808
|
522 |
(*_flow)[e] = 0;
|
kpeter@808
|
523 |
}
|
kpeter@808
|
524 |
}
|
kpeter@808
|
525 |
|
kpeter@808
|
526 |
// Find active nodes (i.e. nodes with positive excess)
|
kpeter@808
|
527 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
kpeter@808
|
528 |
if (_excess[n] > 0) active_nodes.push_back(n);
|
kpeter@808
|
529 |
}
|
kpeter@808
|
530 |
|
kpeter@808
|
531 |
// Initialize the next arc maps
|
kpeter@808
|
532 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
kpeter@808
|
533 |
next_out[n] = OutArcIt(_graph, n);
|
kpeter@808
|
534 |
next_in[n] = InArcIt(_graph, n);
|
kpeter@808
|
535 |
next_dir[n] = true;
|
kpeter@808
|
536 |
}
|
kpeter@808
|
537 |
|
kpeter@808
|
538 |
// Perform partial augment and relabel operations
|
kpeter@808
|
539 |
while (active_nodes.size() > 0) {
|
kpeter@808
|
540 |
// Select an active node (FIFO selection)
|
kpeter@808
|
541 |
if (_excess[active_nodes[0]] <= 0) {
|
kpeter@808
|
542 |
active_nodes.pop_front();
|
kpeter@808
|
543 |
continue;
|
kpeter@808
|
544 |
}
|
kpeter@808
|
545 |
Node start = active_nodes[0];
|
kpeter@808
|
546 |
path_nodes.clear();
|
kpeter@808
|
547 |
path_nodes.push_back(start);
|
kpeter@808
|
548 |
|
kpeter@808
|
549 |
// Find an augmenting path from the start node
|
kpeter@808
|
550 |
Node u, tip = start;
|
kpeter@808
|
551 |
LCost min_red_cost;
|
kpeter@808
|
552 |
while ( _excess[tip] >= 0 &&
|
kpeter@808
|
553 |
int(path_nodes.size()) <= MAX_PATH_LENGTH )
|
kpeter@808
|
554 |
{
|
kpeter@808
|
555 |
if (next_dir[tip]) {
|
kpeter@808
|
556 |
for (OutArcIt e = next_out[tip]; e != INVALID; ++e) {
|
kpeter@808
|
557 |
if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
|
kpeter@808
|
558 |
u = _graph.target(e);
|
kpeter@808
|
559 |
pred_arc[u] = e;
|
kpeter@808
|
560 |
forward[u] = true;
|
kpeter@808
|
561 |
next_out[tip] = e;
|
kpeter@808
|
562 |
tip = u;
|
kpeter@808
|
563 |
path_nodes.push_back(tip);
|
kpeter@808
|
564 |
goto next_step;
|
kpeter@808
|
565 |
}
|
kpeter@808
|
566 |
}
|
kpeter@808
|
567 |
next_dir[tip] = false;
|
kpeter@808
|
568 |
}
|
kpeter@808
|
569 |
for (InArcIt e = next_in[tip]; e != INVALID; ++e) {
|
kpeter@808
|
570 |
if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
|
kpeter@808
|
571 |
u = _graph.source(e);
|
kpeter@808
|
572 |
pred_arc[u] = e;
|
kpeter@808
|
573 |
forward[u] = false;
|
kpeter@808
|
574 |
next_in[tip] = e;
|
kpeter@808
|
575 |
tip = u;
|
kpeter@808
|
576 |
path_nodes.push_back(tip);
|
kpeter@808
|
577 |
goto next_step;
|
kpeter@808
|
578 |
}
|
kpeter@808
|
579 |
}
|
kpeter@808
|
580 |
|
kpeter@808
|
581 |
// Relabel tip node
|
kpeter@808
|
582 |
min_red_cost = std::numeric_limits<LCost>::max() / 2;
|
kpeter@808
|
583 |
for (OutArcIt oe(_graph, tip); oe != INVALID; ++oe) {
|
kpeter@808
|
584 |
if ( _capacity[oe] - (*_flow)[oe] > 0 &&
|
kpeter@808
|
585 |
(*_red_cost)[oe] < min_red_cost )
|
kpeter@808
|
586 |
min_red_cost = (*_red_cost)[oe];
|
kpeter@808
|
587 |
}
|
kpeter@808
|
588 |
for (InArcIt ie(_graph, tip); ie != INVALID; ++ie) {
|
kpeter@808
|
589 |
if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
|
kpeter@808
|
590 |
min_red_cost = -(*_red_cost)[ie];
|
kpeter@808
|
591 |
}
|
kpeter@808
|
592 |
(*_potential)[tip] -= min_red_cost + _epsilon;
|
kpeter@808
|
593 |
|
kpeter@808
|
594 |
// Reset the next arc maps
|
kpeter@808
|
595 |
next_out[tip] = OutArcIt(_graph, tip);
|
kpeter@808
|
596 |
next_in[tip] = InArcIt(_graph, tip);
|
kpeter@808
|
597 |
next_dir[tip] = true;
|
kpeter@808
|
598 |
|
kpeter@808
|
599 |
// Step back
|
kpeter@808
|
600 |
if (tip != start) {
|
kpeter@808
|
601 |
path_nodes.pop_back();
|
kpeter@808
|
602 |
tip = path_nodes[path_nodes.size()-1];
|
kpeter@808
|
603 |
}
|
kpeter@808
|
604 |
|
kpeter@808
|
605 |
next_step:
|
kpeter@808
|
606 |
continue;
|
kpeter@808
|
607 |
}
|
kpeter@808
|
608 |
|
kpeter@808
|
609 |
// Augment along the found path (as much flow as possible)
|
kpeter@808
|
610 |
Capacity delta;
|
kpeter@808
|
611 |
for (int i = 1; i < int(path_nodes.size()); ++i) {
|
kpeter@808
|
612 |
u = path_nodes[i];
|
kpeter@808
|
613 |
delta = forward[u] ?
|
kpeter@808
|
614 |
_capacity[pred_arc[u]] - (*_flow)[pred_arc[u]] :
|
kpeter@808
|
615 |
(*_flow)[pred_arc[u]];
|
kpeter@808
|
616 |
delta = std::min(delta, _excess[path_nodes[i-1]]);
|
kpeter@808
|
617 |
(*_flow)[pred_arc[u]] += forward[u] ? delta : -delta;
|
kpeter@808
|
618 |
_excess[path_nodes[i-1]] -= delta;
|
kpeter@808
|
619 |
_excess[u] += delta;
|
kpeter@808
|
620 |
if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
|
kpeter@808
|
621 |
}
|
kpeter@808
|
622 |
}
|
kpeter@808
|
623 |
}
|
kpeter@808
|
624 |
|
kpeter@808
|
625 |
// Compute node potentials for the original costs
|
kpeter@808
|
626 |
ResidualCostMap<CostMap> res_cost(_orig_cost);
|
kpeter@808
|
627 |
BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
|
kpeter@808
|
628 |
bf(*_res_graph, res_cost);
|
kpeter@808
|
629 |
bf.init(0); bf.start();
|
kpeter@808
|
630 |
for (NodeIt n(_graph); n != INVALID; ++n)
|
kpeter@808
|
631 |
(*_potential)[n] = bf.dist(n);
|
kpeter@808
|
632 |
|
kpeter@808
|
633 |
// Handle non-zero lower bounds
|
kpeter@808
|
634 |
if (_lower) {
|
kpeter@808
|
635 |
for (ArcIt e(_graph); e != INVALID; ++e)
|
kpeter@808
|
636 |
(*_flow)[e] += (*_lower)[e];
|
kpeter@808
|
637 |
}
|
kpeter@808
|
638 |
return true;
|
kpeter@808
|
639 |
}
|
kpeter@808
|
640 |
|
kpeter@808
|
641 |
/// Execute the algorithm performing push and relabel operations.
|
kpeter@808
|
642 |
bool startPushRelabel() {
|
kpeter@808
|
643 |
// Paramters for heuristics
|
kpeter@808
|
644 |
// const int BF_HEURISTIC_EPSILON_BOUND = 1000;
|
kpeter@808
|
645 |
// const int BF_HEURISTIC_BOUND_FACTOR = 3;
|
kpeter@808
|
646 |
|
kpeter@808
|
647 |
typename Digraph::template NodeMap<bool> hyper(_graph, false);
|
kpeter@808
|
648 |
typename Digraph::template NodeMap<Arc> pred_arc(_graph);
|
kpeter@808
|
649 |
typename Digraph::template NodeMap<bool> forward(_graph);
|
kpeter@808
|
650 |
typename Digraph::template NodeMap<OutArcIt> next_out(_graph);
|
kpeter@808
|
651 |
typename Digraph::template NodeMap<InArcIt> next_in(_graph);
|
kpeter@808
|
652 |
typename Digraph::template NodeMap<bool> next_dir(_graph);
|
kpeter@808
|
653 |
std::deque<Node> active_nodes;
|
kpeter@808
|
654 |
|
kpeter@808
|
655 |
// int node_num = countNodes(_graph);
|
kpeter@808
|
656 |
for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
|
kpeter@808
|
657 |
1 : _epsilon / _alpha )
|
kpeter@808
|
658 |
{
|
kpeter@808
|
659 |
/*
|
kpeter@808
|
660 |
// "Early Termination" heuristic: use Bellman-Ford algorithm
|
kpeter@808
|
661 |
// to check if the current flow is optimal
|
kpeter@808
|
662 |
if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
|
kpeter@808
|
663 |
typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
|
kpeter@808
|
664 |
ShiftCostMap shift_cost(_res_cost, 1);
|
kpeter@808
|
665 |
BellmanFord<ResDigraph, ShiftCostMap> bf(*_res_graph, shift_cost);
|
kpeter@808
|
666 |
bf.init(0);
|
kpeter@808
|
667 |
bool done = false;
|
kpeter@808
|
668 |
int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
|
kpeter@808
|
669 |
for (int i = 0; i < K && !done; ++i)
|
kpeter@808
|
670 |
done = bf.processNextWeakRound();
|
kpeter@808
|
671 |
if (done) break;
|
kpeter@808
|
672 |
}
|
kpeter@808
|
673 |
*/
|
kpeter@808
|
674 |
|
kpeter@808
|
675 |
// Saturate arcs not satisfying the optimality condition
|
kpeter@808
|
676 |
Capacity delta;
|
kpeter@808
|
677 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
kpeter@808
|
678 |
if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
|
kpeter@808
|
679 |
delta = _capacity[e] - (*_flow)[e];
|
kpeter@808
|
680 |
_excess[_graph.source(e)] -= delta;
|
kpeter@808
|
681 |
_excess[_graph.target(e)] += delta;
|
kpeter@808
|
682 |
(*_flow)[e] = _capacity[e];
|
kpeter@808
|
683 |
}
|
kpeter@808
|
684 |
if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
|
kpeter@808
|
685 |
_excess[_graph.target(e)] -= (*_flow)[e];
|
kpeter@808
|
686 |
_excess[_graph.source(e)] += (*_flow)[e];
|
kpeter@808
|
687 |
(*_flow)[e] = 0;
|
kpeter@808
|
688 |
}
|
kpeter@808
|
689 |
}
|
kpeter@808
|
690 |
|
kpeter@808
|
691 |
// Find active nodes (i.e. nodes with positive excess)
|
kpeter@808
|
692 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
kpeter@808
|
693 |
if (_excess[n] > 0) active_nodes.push_back(n);
|
kpeter@808
|
694 |
}
|
kpeter@808
|
695 |
|
kpeter@808
|
696 |
// Initialize the next arc maps
|
kpeter@808
|
697 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
kpeter@808
|
698 |
next_out[n] = OutArcIt(_graph, n);
|
kpeter@808
|
699 |
next_in[n] = InArcIt(_graph, n);
|
kpeter@808
|
700 |
next_dir[n] = true;
|
kpeter@808
|
701 |
}
|
kpeter@808
|
702 |
|
kpeter@808
|
703 |
// Perform push and relabel operations
|
kpeter@808
|
704 |
while (active_nodes.size() > 0) {
|
kpeter@808
|
705 |
// Select an active node (FIFO selection)
|
kpeter@808
|
706 |
Node n = active_nodes[0], t;
|
kpeter@808
|
707 |
bool relabel_enabled = true;
|
kpeter@808
|
708 |
|
kpeter@808
|
709 |
// Perform push operations if there are admissible arcs
|
kpeter@808
|
710 |
if (_excess[n] > 0 && next_dir[n]) {
|
kpeter@808
|
711 |
OutArcIt e = next_out[n];
|
kpeter@808
|
712 |
for ( ; e != INVALID; ++e) {
|
kpeter@808
|
713 |
if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
|
kpeter@808
|
714 |
delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
|
kpeter@808
|
715 |
t = _graph.target(e);
|
kpeter@808
|
716 |
|
kpeter@808
|
717 |
// Push-look-ahead heuristic
|
kpeter@808
|
718 |
Capacity ahead = -_excess[t];
|
kpeter@808
|
719 |
for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
|
kpeter@808
|
720 |
if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
|
kpeter@808
|
721 |
ahead += _capacity[oe] - (*_flow)[oe];
|
kpeter@808
|
722 |
}
|
kpeter@808
|
723 |
for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
|
kpeter@808
|
724 |
if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
|
kpeter@808
|
725 |
ahead += (*_flow)[ie];
|
kpeter@808
|
726 |
}
|
kpeter@808
|
727 |
if (ahead < 0) ahead = 0;
|
kpeter@808
|
728 |
|
kpeter@808
|
729 |
// Push flow along the arc
|
kpeter@808
|
730 |
if (ahead < delta) {
|
kpeter@808
|
731 |
(*_flow)[e] += ahead;
|
kpeter@808
|
732 |
_excess[n] -= ahead;
|
kpeter@808
|
733 |
_excess[t] += ahead;
|
kpeter@808
|
734 |
active_nodes.push_front(t);
|
kpeter@808
|
735 |
hyper[t] = true;
|
kpeter@808
|
736 |
relabel_enabled = false;
|
kpeter@808
|
737 |
break;
|
kpeter@808
|
738 |
} else {
|
kpeter@808
|
739 |
(*_flow)[e] += delta;
|
kpeter@808
|
740 |
_excess[n] -= delta;
|
kpeter@808
|
741 |
_excess[t] += delta;
|
kpeter@808
|
742 |
if (_excess[t] > 0 && _excess[t] <= delta)
|
kpeter@808
|
743 |
active_nodes.push_back(t);
|
kpeter@808
|
744 |
}
|
kpeter@808
|
745 |
|
kpeter@808
|
746 |
if (_excess[n] == 0) break;
|
kpeter@808
|
747 |
}
|
kpeter@808
|
748 |
}
|
kpeter@808
|
749 |
if (e != INVALID) {
|
kpeter@808
|
750 |
next_out[n] = e;
|
kpeter@808
|
751 |
} else {
|
kpeter@808
|
752 |
next_dir[n] = false;
|
kpeter@808
|
753 |
}
|
kpeter@808
|
754 |
}
|
kpeter@808
|
755 |
|
kpeter@808
|
756 |
if (_excess[n] > 0 && !next_dir[n]) {
|
kpeter@808
|
757 |
InArcIt e = next_in[n];
|
kpeter@808
|
758 |
for ( ; e != INVALID; ++e) {
|
kpeter@808
|
759 |
if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
|
kpeter@808
|
760 |
delta = std::min((*_flow)[e], _excess[n]);
|
kpeter@808
|
761 |
t = _graph.source(e);
|
kpeter@808
|
762 |
|
kpeter@808
|
763 |
// Push-look-ahead heuristic
|
kpeter@808
|
764 |
Capacity ahead = -_excess[t];
|
kpeter@808
|
765 |
for (OutArcIt oe(_graph, t); oe != INVALID; ++oe) {
|
kpeter@808
|
766 |
if (_capacity[oe] - (*_flow)[oe] > 0 && (*_red_cost)[oe] < 0)
|
kpeter@808
|
767 |
ahead += _capacity[oe] - (*_flow)[oe];
|
kpeter@808
|
768 |
}
|
kpeter@808
|
769 |
for (InArcIt ie(_graph, t); ie != INVALID; ++ie) {
|
kpeter@808
|
770 |
if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < 0)
|
kpeter@808
|
771 |
ahead += (*_flow)[ie];
|
kpeter@808
|
772 |
}
|
kpeter@808
|
773 |
if (ahead < 0) ahead = 0;
|
kpeter@808
|
774 |
|
kpeter@808
|
775 |
// Push flow along the arc
|
kpeter@808
|
776 |
if (ahead < delta) {
|
kpeter@808
|
777 |
(*_flow)[e] -= ahead;
|
kpeter@808
|
778 |
_excess[n] -= ahead;
|
kpeter@808
|
779 |
_excess[t] += ahead;
|
kpeter@808
|
780 |
active_nodes.push_front(t);
|
kpeter@808
|
781 |
hyper[t] = true;
|
kpeter@808
|
782 |
relabel_enabled = false;
|
kpeter@808
|
783 |
break;
|
kpeter@808
|
784 |
} else {
|
kpeter@808
|
785 |
(*_flow)[e] -= delta;
|
kpeter@808
|
786 |
_excess[n] -= delta;
|
kpeter@808
|
787 |
_excess[t] += delta;
|
kpeter@808
|
788 |
if (_excess[t] > 0 && _excess[t] <= delta)
|
kpeter@808
|
789 |
active_nodes.push_back(t);
|
kpeter@808
|
790 |
}
|
kpeter@808
|
791 |
|
kpeter@808
|
792 |
if (_excess[n] == 0) break;
|
kpeter@808
|
793 |
}
|
kpeter@808
|
794 |
}
|
kpeter@808
|
795 |
next_in[n] = e;
|
kpeter@808
|
796 |
}
|
kpeter@808
|
797 |
|
kpeter@808
|
798 |
// Relabel the node if it is still active (or hyper)
|
kpeter@808
|
799 |
if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
|
kpeter@808
|
800 |
LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
|
kpeter@808
|
801 |
for (OutArcIt oe(_graph, n); oe != INVALID; ++oe) {
|
kpeter@808
|
802 |
if ( _capacity[oe] - (*_flow)[oe] > 0 &&
|
kpeter@808
|
803 |
(*_red_cost)[oe] < min_red_cost )
|
kpeter@808
|
804 |
min_red_cost = (*_red_cost)[oe];
|
kpeter@808
|
805 |
}
|
kpeter@808
|
806 |
for (InArcIt ie(_graph, n); ie != INVALID; ++ie) {
|
kpeter@808
|
807 |
if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
|
kpeter@808
|
808 |
min_red_cost = -(*_red_cost)[ie];
|
kpeter@808
|
809 |
}
|
kpeter@808
|
810 |
(*_potential)[n] -= min_red_cost + _epsilon;
|
kpeter@808
|
811 |
hyper[n] = false;
|
kpeter@808
|
812 |
|
kpeter@808
|
813 |
// Reset the next arc maps
|
kpeter@808
|
814 |
next_out[n] = OutArcIt(_graph, n);
|
kpeter@808
|
815 |
next_in[n] = InArcIt(_graph, n);
|
kpeter@808
|
816 |
next_dir[n] = true;
|
kpeter@808
|
817 |
}
|
kpeter@808
|
818 |
|
kpeter@808
|
819 |
// Remove nodes that are not active nor hyper
|
kpeter@808
|
820 |
while ( active_nodes.size() > 0 &&
|
kpeter@808
|
821 |
_excess[active_nodes[0]] <= 0 &&
|
kpeter@808
|
822 |
!hyper[active_nodes[0]] ) {
|
kpeter@808
|
823 |
active_nodes.pop_front();
|
kpeter@808
|
824 |
}
|
kpeter@808
|
825 |
}
|
kpeter@808
|
826 |
}
|
kpeter@808
|
827 |
|
kpeter@808
|
828 |
// Compute node potentials for the original costs
|
kpeter@808
|
829 |
ResidualCostMap<CostMap> res_cost(_orig_cost);
|
kpeter@808
|
830 |
BellmanFord< ResDigraph, ResidualCostMap<CostMap> >
|
kpeter@808
|
831 |
bf(*_res_graph, res_cost);
|
kpeter@808
|
832 |
bf.init(0); bf.start();
|
kpeter@808
|
833 |
for (NodeIt n(_graph); n != INVALID; ++n)
|
kpeter@808
|
834 |
(*_potential)[n] = bf.dist(n);
|
kpeter@808
|
835 |
|
kpeter@808
|
836 |
// Handle non-zero lower bounds
|
kpeter@808
|
837 |
if (_lower) {
|
kpeter@808
|
838 |
for (ArcIt e(_graph); e != INVALID; ++e)
|
kpeter@808
|
839 |
(*_flow)[e] += (*_lower)[e];
|
kpeter@808
|
840 |
}
|
kpeter@808
|
841 |
return true;
|
kpeter@808
|
842 |
}
|
kpeter@808
|
843 |
|
kpeter@808
|
844 |
}; //class CostScaling
|
kpeter@808
|
845 |
|
kpeter@808
|
846 |
///@}
|
kpeter@808
|
847 |
|
kpeter@808
|
848 |
} //namespace lemon
|
kpeter@808
|
849 |
|
kpeter@808
|
850 |
#endif //LEMON_COST_SCALING_H
|