1 |
1 |
/* -*- mode: C++; indent-tabs-mode: nil; -*-
|
2 |
2 |
*
|
3 |
3 |
* This file is a part of LEMON, a generic C++ optimization library.
|
4 |
4 |
*
|
5 |
5 |
* Copyright (C) 2003-2009
|
6 |
6 |
* Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
|
7 |
7 |
* (Egervary Research Group on Combinatorial Optimization, EGRES).
|
8 |
8 |
*
|
9 |
9 |
* Permission to use, modify and distribute this software is granted
|
10 |
10 |
* provided that this copyright notice appears in all copies. For
|
11 |
11 |
* precise terms see the accompanying LICENSE file.
|
12 |
12 |
*
|
13 |
13 |
* This software is provided "AS IS" with no warranty of any kind,
|
14 |
14 |
* express or implied, and with no claim as to its suitability for any
|
15 |
15 |
* purpose.
|
16 |
16 |
*
|
17 |
17 |
*/
|
18 |
18 |
|
19 |
19 |
#ifndef LEMON_NETWORK_SIMPLEX_H
|
20 |
20 |
#define LEMON_NETWORK_SIMPLEX_H
|
21 |
21 |
|
22 |
22 |
/// \ingroup min_cost_flow
|
23 |
23 |
///
|
24 |
24 |
/// \file
|
25 |
25 |
/// \brief Network Simplex algorithm for finding a minimum cost flow.
|
26 |
26 |
|
27 |
27 |
#include <vector>
|
28 |
28 |
#include <limits>
|
29 |
29 |
#include <algorithm>
|
30 |
30 |
|
31 |
31 |
#include <lemon/core.h>
|
32 |
32 |
#include <lemon/math.h>
|
|
33 |
#include <lemon/maps.h>
|
|
34 |
#include <lemon/circulation.h>
|
|
35 |
#include <lemon/adaptors.h>
|
33 |
36 |
|
34 |
37 |
namespace lemon {
|
35 |
38 |
|
36 |
39 |
/// \addtogroup min_cost_flow
|
37 |
40 |
/// @{
|
38 |
41 |
|
39 |
42 |
/// \brief Implementation of the primal Network Simplex algorithm
|
40 |
43 |
/// for finding a \ref min_cost_flow "minimum cost flow".
|
41 |
44 |
///
|
42 |
45 |
/// \ref NetworkSimplex implements the primal Network Simplex algorithm
|
43 |
46 |
/// for finding a \ref min_cost_flow "minimum cost flow".
|
44 |
47 |
/// This algorithm is a specialized version of the linear programming
|
45 |
48 |
/// simplex method directly for the minimum cost flow problem.
|
46 |
49 |
/// It is one of the most efficient solution methods.
|
47 |
50 |
///
|
48 |
51 |
/// In general this class is the fastest implementation available
|
49 |
52 |
/// in LEMON for the minimum cost flow problem.
|
|
53 |
/// Moreover it supports both direction of the supply/demand inequality
|
|
54 |
/// constraints. For more information see \ref ProblemType.
|
50 |
55 |
///
|
51 |
56 |
/// \tparam GR The digraph type the algorithm runs on.
|
52 |
57 |
/// \tparam F The value type used for flow amounts, capacity bounds
|
53 |
58 |
/// and supply values in the algorithm. By default it is \c int.
|
54 |
59 |
/// \tparam C The value type used for costs and potentials in the
|
55 |
60 |
/// algorithm. By default it is the same as \c F.
|
56 |
61 |
///
|
57 |
62 |
/// \warning Both value types must be signed and all input data must
|
58 |
63 |
/// be integer.
|
59 |
64 |
///
|
60 |
65 |
/// \note %NetworkSimplex provides five different pivot rule
|
61 |
|
/// implementations. For more information see \ref PivotRule.
|
|
66 |
/// implementations, from which the most efficient one is used
|
|
67 |
/// by default. For more information see \ref PivotRule.
|
62 |
68 |
template <typename GR, typename F = int, typename C = F>
|
63 |
69 |
class NetworkSimplex
|
64 |
70 |
{
|
65 |
71 |
public:
|
66 |
72 |
|
67 |
73 |
/// The flow type of the algorithm
|
68 |
74 |
typedef F Flow;
|
69 |
75 |
/// The cost type of the algorithm
|
70 |
76 |
typedef C Cost;
|
|
77 |
#ifdef DOXYGEN
|
|
78 |
/// The type of the flow map
|
|
79 |
typedef GR::ArcMap<Flow> FlowMap;
|
|
80 |
/// The type of the potential map
|
|
81 |
typedef GR::NodeMap<Cost> PotentialMap;
|
|
82 |
#else
|
71 |
83 |
/// The type of the flow map
|
72 |
84 |
typedef typename GR::template ArcMap<Flow> FlowMap;
|
73 |
85 |
/// The type of the potential map
|
74 |
86 |
typedef typename GR::template NodeMap<Cost> PotentialMap;
|
|
87 |
#endif
|
75 |
88 |
|
76 |
89 |
public:
|
77 |
90 |
|
78 |
91 |
/// \brief Enum type for selecting the pivot rule.
|
79 |
92 |
///
|
80 |
93 |
/// Enum type for selecting the pivot rule for the \ref run()
|
81 |
94 |
/// function.
|
82 |
95 |
///
|
83 |
96 |
/// \ref NetworkSimplex provides five different pivot rule
|
84 |
97 |
/// implementations that significantly affect the running time
|
85 |
98 |
/// of the algorithm.
|
86 |
99 |
/// By default \ref BLOCK_SEARCH "Block Search" is used, which
|
87 |
100 |
/// proved to be the most efficient and the most robust on various
|
88 |
101 |
/// test inputs according to our benchmark tests.
|
89 |
102 |
/// However another pivot rule can be selected using the \ref run()
|
90 |
103 |
/// function with the proper parameter.
|
91 |
104 |
enum PivotRule {
|
92 |
105 |
|
93 |
106 |
/// The First Eligible pivot rule.
|
94 |
107 |
/// The next eligible arc is selected in a wraparound fashion
|
95 |
108 |
/// in every iteration.
|
96 |
109 |
FIRST_ELIGIBLE,
|
97 |
110 |
|
98 |
111 |
/// The Best Eligible pivot rule.
|
99 |
112 |
/// The best eligible arc is selected in every iteration.
|
100 |
113 |
BEST_ELIGIBLE,
|
101 |
114 |
|
102 |
115 |
/// The Block Search pivot rule.
|
103 |
116 |
/// A specified number of arcs are examined in every iteration
|
104 |
117 |
/// in a wraparound fashion and the best eligible arc is selected
|
105 |
118 |
/// from this block.
|
106 |
119 |
BLOCK_SEARCH,
|
107 |
120 |
|
108 |
121 |
/// The Candidate List pivot rule.
|
109 |
122 |
/// In a major iteration a candidate list is built from eligible arcs
|
110 |
123 |
/// in a wraparound fashion and in the following minor iterations
|
111 |
124 |
/// the best eligible arc is selected from this list.
|
112 |
125 |
CANDIDATE_LIST,
|
113 |
126 |
|
114 |
127 |
/// The Altering Candidate List pivot rule.
|
115 |
128 |
/// It is a modified version of the Candidate List method.
|
116 |
129 |
/// It keeps only the several best eligible arcs from the former
|
117 |
130 |
/// candidate list and extends this list in every iteration.
|
118 |
131 |
ALTERING_LIST
|
119 |
132 |
};
|
|
133 |
|
|
134 |
/// \brief Enum type for selecting the problem type.
|
|
135 |
///
|
|
136 |
/// Enum type for selecting the problem type, i.e. the direction of
|
|
137 |
/// the inequalities in the supply/demand constraints of the
|
|
138 |
/// \ref min_cost_flow "minimum cost flow problem".
|
|
139 |
///
|
|
140 |
/// The default problem type is \c GEQ, since this form is supported
|
|
141 |
/// by other minimum cost flow algorithms and the \ref Circulation
|
|
142 |
/// algorithm as well.
|
|
143 |
/// The \c LEQ problem type can be selected using the \ref problemType()
|
|
144 |
/// function.
|
|
145 |
///
|
|
146 |
/// Note that the equality form is a special case of both problem type.
|
|
147 |
enum ProblemType {
|
|
148 |
|
|
149 |
/// This option means that there are "<em>greater or equal</em>"
|
|
150 |
/// constraints in the defintion, i.e. the exact formulation of the
|
|
151 |
/// problem is the following.
|
|
152 |
/**
|
|
153 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
|
|
154 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
|
|
155 |
sup(u) \quad \forall u\in V \f]
|
|
156 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
|
|
157 |
*/
|
|
158 |
/// It means that the total demand must be greater or equal to the
|
|
159 |
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
|
|
160 |
/// negative) and all the supplies have to be carried out from
|
|
161 |
/// the supply nodes, but there could be demands that are not
|
|
162 |
/// satisfied.
|
|
163 |
GEQ,
|
|
164 |
/// It is just an alias for the \c GEQ option.
|
|
165 |
CARRY_SUPPLIES = GEQ,
|
|
166 |
|
|
167 |
/// This option means that there are "<em>less or equal</em>"
|
|
168 |
/// constraints in the defintion, i.e. the exact formulation of the
|
|
169 |
/// problem is the following.
|
|
170 |
/**
|
|
171 |
\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
|
|
172 |
\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
|
|
173 |
sup(u) \quad \forall u\in V \f]
|
|
174 |
\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
|
|
175 |
*/
|
|
176 |
/// It means that the total demand must be less or equal to the
|
|
177 |
/// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
|
|
178 |
/// positive) and all the demands have to be satisfied, but there
|
|
179 |
/// could be supplies that are not carried out from the supply
|
|
180 |
/// nodes.
|
|
181 |
LEQ,
|
|
182 |
/// It is just an alias for the \c LEQ option.
|
|
183 |
SATISFY_DEMANDS = LEQ
|
|
184 |
};
|
120 |
185 |
|
121 |
186 |
private:
|
122 |
187 |
|
123 |
188 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR);
|
124 |
189 |
|
125 |
190 |
typedef typename GR::template ArcMap<Flow> FlowArcMap;
|
126 |
191 |
typedef typename GR::template ArcMap<Cost> CostArcMap;
|
127 |
192 |
typedef typename GR::template NodeMap<Flow> FlowNodeMap;
|
128 |
193 |
|
129 |
194 |
typedef std::vector<Arc> ArcVector;
|
130 |
195 |
typedef std::vector<Node> NodeVector;
|
131 |
196 |
typedef std::vector<int> IntVector;
|
132 |
197 |
typedef std::vector<bool> BoolVector;
|
133 |
198 |
typedef std::vector<Flow> FlowVector;
|
134 |
199 |
typedef std::vector<Cost> CostVector;
|
135 |
200 |
|
136 |
201 |
// State constants for arcs
|
137 |
202 |
enum ArcStateEnum {
|
138 |
203 |
STATE_UPPER = -1,
|
139 |
204 |
STATE_TREE = 0,
|
140 |
205 |
STATE_LOWER = 1
|
141 |
206 |
};
|
142 |
207 |
|
143 |
208 |
private:
|
144 |
209 |
|
145 |
210 |
// Data related to the underlying digraph
|
146 |
211 |
const GR &_graph;
|
147 |
212 |
int _node_num;
|
148 |
213 |
int _arc_num;
|
149 |
214 |
|
150 |
215 |
// Parameters of the problem
|
151 |
216 |
FlowArcMap *_plower;
|
152 |
217 |
FlowArcMap *_pupper;
|
153 |
218 |
CostArcMap *_pcost;
|
154 |
219 |
FlowNodeMap *_psupply;
|
155 |
220 |
bool _pstsup;
|
156 |
221 |
Node _psource, _ptarget;
|
157 |
222 |
Flow _pstflow;
|
|
223 |
ProblemType _ptype;
|
158 |
224 |
|
159 |
225 |
// Result maps
|
160 |
226 |
FlowMap *_flow_map;
|
161 |
227 |
PotentialMap *_potential_map;
|
162 |
228 |
bool _local_flow;
|
163 |
229 |
bool _local_potential;
|
164 |
230 |
|
165 |
231 |
// Data structures for storing the digraph
|
166 |
232 |
IntNodeMap _node_id;
|
167 |
233 |
ArcVector _arc_ref;
|
168 |
234 |
IntVector _source;
|
169 |
235 |
IntVector _target;
|
170 |
236 |
|
171 |
237 |
// Node and arc data
|
172 |
238 |
FlowVector _cap;
|
173 |
239 |
CostVector _cost;
|
174 |
240 |
FlowVector _supply;
|
175 |
241 |
FlowVector _flow;
|
176 |
242 |
CostVector _pi;
|
177 |
243 |
|
178 |
244 |
// Data for storing the spanning tree structure
|
179 |
245 |
IntVector _parent;
|
180 |
246 |
IntVector _pred;
|
181 |
247 |
IntVector _thread;
|
182 |
248 |
IntVector _rev_thread;
|
183 |
249 |
IntVector _succ_num;
|
184 |
250 |
IntVector _last_succ;
|
185 |
251 |
IntVector _dirty_revs;
|
186 |
252 |
BoolVector _forward;
|
187 |
253 |
IntVector _state;
|
188 |
254 |
int _root;
|
189 |
255 |
|
190 |
256 |
// Temporary data used in the current pivot iteration
|
191 |
257 |
int in_arc, join, u_in, v_in, u_out, v_out;
|
192 |
258 |
int first, second, right, last;
|
193 |
259 |
int stem, par_stem, new_stem;
|
194 |
260 |
Flow delta;
|
195 |
261 |
|
196 |
262 |
private:
|
197 |
263 |
|
198 |
264 |
// Implementation of the First Eligible pivot rule
|
199 |
265 |
class FirstEligiblePivotRule
|
200 |
266 |
{
|
201 |
267 |
private:
|
202 |
268 |
|
203 |
269 |
// References to the NetworkSimplex class
|
204 |
270 |
const IntVector &_source;
|
205 |
271 |
const IntVector &_target;
|
206 |
272 |
const CostVector &_cost;
|
207 |
273 |
const IntVector &_state;
|
208 |
274 |
const CostVector &_pi;
|
209 |
275 |
int &_in_arc;
|
210 |
276 |
int _arc_num;
|
211 |
277 |
|
212 |
278 |
// Pivot rule data
|
213 |
279 |
int _next_arc;
|
214 |
280 |
|
215 |
281 |
public:
|
216 |
282 |
|
217 |
283 |
// Constructor
|
218 |
284 |
FirstEligiblePivotRule(NetworkSimplex &ns) :
|
219 |
285 |
_source(ns._source), _target(ns._target),
|
220 |
286 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
221 |
287 |
_in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
|
222 |
288 |
{}
|
223 |
289 |
|
224 |
290 |
// Find next entering arc
|
225 |
291 |
bool findEnteringArc() {
|
226 |
292 |
Cost c;
|
227 |
293 |
for (int e = _next_arc; e < _arc_num; ++e) {
|
228 |
294 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
229 |
295 |
if (c < 0) {
|
230 |
296 |
_in_arc = e;
|
231 |
297 |
_next_arc = e + 1;
|
232 |
298 |
return true;
|
233 |
299 |
}
|
234 |
300 |
}
|
235 |
301 |
for (int e = 0; e < _next_arc; ++e) {
|
236 |
302 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
237 |
303 |
if (c < 0) {
|
238 |
304 |
_in_arc = e;
|
239 |
305 |
_next_arc = e + 1;
|
240 |
306 |
return true;
|
241 |
307 |
}
|
242 |
308 |
}
|
243 |
309 |
return false;
|
244 |
310 |
}
|
245 |
311 |
|
246 |
312 |
}; //class FirstEligiblePivotRule
|
247 |
313 |
|
248 |
314 |
|
249 |
315 |
// Implementation of the Best Eligible pivot rule
|
250 |
316 |
class BestEligiblePivotRule
|
251 |
317 |
{
|
252 |
318 |
private:
|
253 |
319 |
|
... |
... |
@@ -493,732 +559,838 @@
|
493 |
559 |
return _map[left] > _map[right];
|
494 |
560 |
}
|
495 |
561 |
};
|
496 |
562 |
|
497 |
563 |
SortFunc _sort_func;
|
498 |
564 |
|
499 |
565 |
public:
|
500 |
566 |
|
501 |
567 |
// Constructor
|
502 |
568 |
AlteringListPivotRule(NetworkSimplex &ns) :
|
503 |
569 |
_source(ns._source), _target(ns._target),
|
504 |
570 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
|
505 |
571 |
_in_arc(ns.in_arc), _arc_num(ns._arc_num),
|
506 |
572 |
_next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
|
507 |
573 |
{
|
508 |
574 |
// The main parameters of the pivot rule
|
509 |
575 |
const double BLOCK_SIZE_FACTOR = 1.5;
|
510 |
576 |
const int MIN_BLOCK_SIZE = 10;
|
511 |
577 |
const double HEAD_LENGTH_FACTOR = 0.1;
|
512 |
578 |
const int MIN_HEAD_LENGTH = 3;
|
513 |
579 |
|
514 |
580 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_arc_num)),
|
515 |
581 |
MIN_BLOCK_SIZE );
|
516 |
582 |
_head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
|
517 |
583 |
MIN_HEAD_LENGTH );
|
518 |
584 |
_candidates.resize(_head_length + _block_size);
|
519 |
585 |
_curr_length = 0;
|
520 |
586 |
}
|
521 |
587 |
|
522 |
588 |
// Find next entering arc
|
523 |
589 |
bool findEnteringArc() {
|
524 |
590 |
// Check the current candidate list
|
525 |
591 |
int e;
|
526 |
592 |
for (int i = 0; i < _curr_length; ++i) {
|
527 |
593 |
e = _candidates[i];
|
528 |
594 |
_cand_cost[e] = _state[e] *
|
529 |
595 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
530 |
596 |
if (_cand_cost[e] >= 0) {
|
531 |
597 |
_candidates[i--] = _candidates[--_curr_length];
|
532 |
598 |
}
|
533 |
599 |
}
|
534 |
600 |
|
535 |
601 |
// Extend the list
|
536 |
602 |
int cnt = _block_size;
|
537 |
603 |
int last_arc = 0;
|
538 |
604 |
int limit = _head_length;
|
539 |
605 |
|
540 |
606 |
for (int e = _next_arc; e < _arc_num; ++e) {
|
541 |
607 |
_cand_cost[e] = _state[e] *
|
542 |
608 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
543 |
609 |
if (_cand_cost[e] < 0) {
|
544 |
610 |
_candidates[_curr_length++] = e;
|
545 |
611 |
last_arc = e;
|
546 |
612 |
}
|
547 |
613 |
if (--cnt == 0) {
|
548 |
614 |
if (_curr_length > limit) break;
|
549 |
615 |
limit = 0;
|
550 |
616 |
cnt = _block_size;
|
551 |
617 |
}
|
552 |
618 |
}
|
553 |
619 |
if (_curr_length <= limit) {
|
554 |
620 |
for (int e = 0; e < _next_arc; ++e) {
|
555 |
621 |
_cand_cost[e] = _state[e] *
|
556 |
622 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
|
557 |
623 |
if (_cand_cost[e] < 0) {
|
558 |
624 |
_candidates[_curr_length++] = e;
|
559 |
625 |
last_arc = e;
|
560 |
626 |
}
|
561 |
627 |
if (--cnt == 0) {
|
562 |
628 |
if (_curr_length > limit) break;
|
563 |
629 |
limit = 0;
|
564 |
630 |
cnt = _block_size;
|
565 |
631 |
}
|
566 |
632 |
}
|
567 |
633 |
}
|
568 |
634 |
if (_curr_length == 0) return false;
|
569 |
635 |
_next_arc = last_arc + 1;
|
570 |
636 |
|
571 |
637 |
// Make heap of the candidate list (approximating a partial sort)
|
572 |
638 |
make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
|
573 |
639 |
_sort_func );
|
574 |
640 |
|
575 |
641 |
// Pop the first element of the heap
|
576 |
642 |
_in_arc = _candidates[0];
|
577 |
643 |
pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
|
578 |
644 |
_sort_func );
|
579 |
645 |
_curr_length = std::min(_head_length, _curr_length - 1);
|
580 |
646 |
return true;
|
581 |
647 |
}
|
582 |
648 |
|
583 |
649 |
}; //class AlteringListPivotRule
|
584 |
650 |
|
585 |
651 |
public:
|
586 |
652 |
|
587 |
653 |
/// \brief Constructor.
|
588 |
654 |
///
|
589 |
|
/// Constructor.
|
|
655 |
/// The constructor of the class.
|
590 |
656 |
///
|
591 |
657 |
/// \param graph The digraph the algorithm runs on.
|
592 |
658 |
NetworkSimplex(const GR& graph) :
|
593 |
659 |
_graph(graph),
|
594 |
660 |
_plower(NULL), _pupper(NULL), _pcost(NULL),
|
595 |
|
_psupply(NULL), _pstsup(false),
|
|
661 |
_psupply(NULL), _pstsup(false), _ptype(GEQ),
|
596 |
662 |
_flow_map(NULL), _potential_map(NULL),
|
597 |
663 |
_local_flow(false), _local_potential(false),
|
598 |
664 |
_node_id(graph)
|
599 |
665 |
{
|
600 |
666 |
LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
|
601 |
667 |
std::numeric_limits<Flow>::is_signed,
|
602 |
668 |
"The flow type of NetworkSimplex must be signed integer");
|
603 |
669 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_integer &&
|
604 |
670 |
std::numeric_limits<Cost>::is_signed,
|
605 |
671 |
"The cost type of NetworkSimplex must be signed integer");
|
606 |
672 |
}
|
607 |
673 |
|
608 |
674 |
/// Destructor.
|
609 |
675 |
~NetworkSimplex() {
|
610 |
676 |
if (_local_flow) delete _flow_map;
|
611 |
677 |
if (_local_potential) delete _potential_map;
|
612 |
678 |
}
|
613 |
679 |
|
|
680 |
/// \name Parameters
|
|
681 |
/// The parameters of the algorithm can be specified using these
|
|
682 |
/// functions.
|
|
683 |
|
|
684 |
/// @{
|
|
685 |
|
614 |
686 |
/// \brief Set the lower bounds on the arcs.
|
615 |
687 |
///
|
616 |
688 |
/// This function sets the lower bounds on the arcs.
|
617 |
689 |
/// If neither this function nor \ref boundMaps() is used before
|
618 |
690 |
/// calling \ref run(), the lower bounds will be set to zero
|
619 |
691 |
/// on all arcs.
|
620 |
692 |
///
|
621 |
693 |
/// \param map An arc map storing the lower bounds.
|
622 |
694 |
/// Its \c Value type must be convertible to the \c Flow type
|
623 |
695 |
/// of the algorithm.
|
624 |
696 |
///
|
625 |
697 |
/// \return <tt>(*this)</tt>
|
626 |
698 |
template <typename LOWER>
|
627 |
699 |
NetworkSimplex& lowerMap(const LOWER& map) {
|
628 |
700 |
delete _plower;
|
629 |
701 |
_plower = new FlowArcMap(_graph);
|
630 |
702 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
631 |
703 |
(*_plower)[a] = map[a];
|
632 |
704 |
}
|
633 |
705 |
return *this;
|
634 |
706 |
}
|
635 |
707 |
|
636 |
708 |
/// \brief Set the upper bounds (capacities) on the arcs.
|
637 |
709 |
///
|
638 |
710 |
/// This function sets the upper bounds (capacities) on the arcs.
|
639 |
711 |
/// If none of the functions \ref upperMap(), \ref capacityMap()
|
640 |
712 |
/// and \ref boundMaps() is used before calling \ref run(),
|
641 |
713 |
/// the upper bounds (capacities) will be set to
|
642 |
714 |
/// \c std::numeric_limits<Flow>::max() on all arcs.
|
643 |
715 |
///
|
644 |
716 |
/// \param map An arc map storing the upper bounds.
|
645 |
717 |
/// Its \c Value type must be convertible to the \c Flow type
|
646 |
718 |
/// of the algorithm.
|
647 |
719 |
///
|
648 |
720 |
/// \return <tt>(*this)</tt>
|
649 |
721 |
template<typename UPPER>
|
650 |
722 |
NetworkSimplex& upperMap(const UPPER& map) {
|
651 |
723 |
delete _pupper;
|
652 |
724 |
_pupper = new FlowArcMap(_graph);
|
653 |
725 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
654 |
726 |
(*_pupper)[a] = map[a];
|
655 |
727 |
}
|
656 |
728 |
return *this;
|
657 |
729 |
}
|
658 |
730 |
|
659 |
731 |
/// \brief Set the upper bounds (capacities) on the arcs.
|
660 |
732 |
///
|
661 |
733 |
/// This function sets the upper bounds (capacities) on the arcs.
|
662 |
734 |
/// It is just an alias for \ref upperMap().
|
663 |
735 |
///
|
664 |
736 |
/// \return <tt>(*this)</tt>
|
665 |
737 |
template<typename CAP>
|
666 |
738 |
NetworkSimplex& capacityMap(const CAP& map) {
|
667 |
739 |
return upperMap(map);
|
668 |
740 |
}
|
669 |
741 |
|
670 |
742 |
/// \brief Set the lower and upper bounds on the arcs.
|
671 |
743 |
///
|
672 |
744 |
/// This function sets the lower and upper bounds on the arcs.
|
673 |
745 |
/// If neither this function nor \ref lowerMap() is used before
|
674 |
746 |
/// calling \ref run(), the lower bounds will be set to zero
|
675 |
747 |
/// on all arcs.
|
676 |
748 |
/// If none of the functions \ref upperMap(), \ref capacityMap()
|
677 |
749 |
/// and \ref boundMaps() is used before calling \ref run(),
|
678 |
750 |
/// the upper bounds (capacities) will be set to
|
679 |
751 |
/// \c std::numeric_limits<Flow>::max() on all arcs.
|
680 |
752 |
///
|
681 |
753 |
/// \param lower An arc map storing the lower bounds.
|
682 |
754 |
/// \param upper An arc map storing the upper bounds.
|
683 |
755 |
///
|
684 |
756 |
/// The \c Value type of the maps must be convertible to the
|
685 |
757 |
/// \c Flow type of the algorithm.
|
686 |
758 |
///
|
687 |
759 |
/// \note This function is just a shortcut of calling \ref lowerMap()
|
688 |
760 |
/// and \ref upperMap() separately.
|
689 |
761 |
///
|
690 |
762 |
/// \return <tt>(*this)</tt>
|
691 |
763 |
template <typename LOWER, typename UPPER>
|
692 |
764 |
NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {
|
693 |
765 |
return lowerMap(lower).upperMap(upper);
|
694 |
766 |
}
|
695 |
767 |
|
696 |
768 |
/// \brief Set the costs of the arcs.
|
697 |
769 |
///
|
698 |
770 |
/// This function sets the costs of the arcs.
|
699 |
771 |
/// If it is not used before calling \ref run(), the costs
|
700 |
772 |
/// will be set to \c 1 on all arcs.
|
701 |
773 |
///
|
702 |
774 |
/// \param map An arc map storing the costs.
|
703 |
775 |
/// Its \c Value type must be convertible to the \c Cost type
|
704 |
776 |
/// of the algorithm.
|
705 |
777 |
///
|
706 |
778 |
/// \return <tt>(*this)</tt>
|
707 |
779 |
template<typename COST>
|
708 |
780 |
NetworkSimplex& costMap(const COST& map) {
|
709 |
781 |
delete _pcost;
|
710 |
782 |
_pcost = new CostArcMap(_graph);
|
711 |
783 |
for (ArcIt a(_graph); a != INVALID; ++a) {
|
712 |
784 |
(*_pcost)[a] = map[a];
|
713 |
785 |
}
|
714 |
786 |
return *this;
|
715 |
787 |
}
|
716 |
788 |
|
717 |
789 |
/// \brief Set the supply values of the nodes.
|
718 |
790 |
///
|
719 |
791 |
/// This function sets the supply values of the nodes.
|
720 |
792 |
/// If neither this function nor \ref stSupply() is used before
|
721 |
793 |
/// calling \ref run(), the supply of each node will be set to zero.
|
722 |
794 |
/// (It makes sense only if non-zero lower bounds are given.)
|
723 |
795 |
///
|
724 |
796 |
/// \param map A node map storing the supply values.
|
725 |
797 |
/// Its \c Value type must be convertible to the \c Flow type
|
726 |
798 |
/// of the algorithm.
|
727 |
799 |
///
|
728 |
800 |
/// \return <tt>(*this)</tt>
|
729 |
801 |
template<typename SUP>
|
730 |
802 |
NetworkSimplex& supplyMap(const SUP& map) {
|
731 |
803 |
delete _psupply;
|
732 |
804 |
_pstsup = false;
|
733 |
805 |
_psupply = new FlowNodeMap(_graph);
|
734 |
806 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
735 |
807 |
(*_psupply)[n] = map[n];
|
736 |
808 |
}
|
737 |
809 |
return *this;
|
738 |
810 |
}
|
739 |
811 |
|
740 |
812 |
/// \brief Set single source and target nodes and a supply value.
|
741 |
813 |
///
|
742 |
814 |
/// This function sets a single source node and a single target node
|
743 |
815 |
/// and the required flow value.
|
744 |
816 |
/// If neither this function nor \ref supplyMap() is used before
|
745 |
817 |
/// calling \ref run(), the supply of each node will be set to zero.
|
746 |
818 |
/// (It makes sense only if non-zero lower bounds are given.)
|
747 |
819 |
///
|
748 |
820 |
/// \param s The source node.
|
749 |
821 |
/// \param t The target node.
|
750 |
822 |
/// \param k The required amount of flow from node \c s to node \c t
|
751 |
823 |
/// (i.e. the supply of \c s and the demand of \c t).
|
752 |
824 |
///
|
753 |
825 |
/// \return <tt>(*this)</tt>
|
754 |
826 |
NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) {
|
755 |
827 |
delete _psupply;
|
756 |
828 |
_psupply = NULL;
|
757 |
829 |
_pstsup = true;
|
758 |
830 |
_psource = s;
|
759 |
831 |
_ptarget = t;
|
760 |
832 |
_pstflow = k;
|
761 |
833 |
return *this;
|
762 |
834 |
}
|
|
835 |
|
|
836 |
/// \brief Set the problem type.
|
|
837 |
///
|
|
838 |
/// This function sets the problem type for the algorithm.
|
|
839 |
/// If it is not used before calling \ref run(), the \ref GEQ problem
|
|
840 |
/// type will be used.
|
|
841 |
///
|
|
842 |
/// For more information see \ref ProblemType.
|
|
843 |
///
|
|
844 |
/// \return <tt>(*this)</tt>
|
|
845 |
NetworkSimplex& problemType(ProblemType problem_type) {
|
|
846 |
_ptype = problem_type;
|
|
847 |
return *this;
|
|
848 |
}
|
763 |
849 |
|
764 |
850 |
/// \brief Set the flow map.
|
765 |
851 |
///
|
766 |
852 |
/// This function sets the flow map.
|
767 |
853 |
/// If it is not used before calling \ref run(), an instance will
|
768 |
854 |
/// be allocated automatically. The destructor deallocates this
|
769 |
855 |
/// automatically allocated map, of course.
|
770 |
856 |
///
|
771 |
857 |
/// \return <tt>(*this)</tt>
|
772 |
858 |
NetworkSimplex& flowMap(FlowMap& map) {
|
773 |
859 |
if (_local_flow) {
|
774 |
860 |
delete _flow_map;
|
775 |
861 |
_local_flow = false;
|
776 |
862 |
}
|
777 |
863 |
_flow_map = ↦
|
778 |
864 |
return *this;
|
779 |
865 |
}
|
780 |
866 |
|
781 |
867 |
/// \brief Set the potential map.
|
782 |
868 |
///
|
783 |
869 |
/// This function sets the potential map, which is used for storing
|
784 |
870 |
/// the dual solution.
|
785 |
871 |
/// If it is not used before calling \ref run(), an instance will
|
786 |
872 |
/// be allocated automatically. The destructor deallocates this
|
787 |
873 |
/// automatically allocated map, of course.
|
788 |
874 |
///
|
789 |
875 |
/// \return <tt>(*this)</tt>
|
790 |
876 |
NetworkSimplex& potentialMap(PotentialMap& map) {
|
791 |
877 |
if (_local_potential) {
|
792 |
878 |
delete _potential_map;
|
793 |
879 |
_local_potential = false;
|
794 |
880 |
}
|
795 |
881 |
_potential_map = ↦
|
796 |
882 |
return *this;
|
797 |
883 |
}
|
|
884 |
|
|
885 |
/// @}
|
798 |
886 |
|
799 |
887 |
/// \name Execution Control
|
800 |
888 |
/// The algorithm can be executed using \ref run().
|
801 |
889 |
|
802 |
890 |
/// @{
|
803 |
891 |
|
804 |
892 |
/// \brief Run the algorithm.
|
805 |
893 |
///
|
806 |
894 |
/// This function runs the algorithm.
|
807 |
|
/// The paramters can be specified using \ref lowerMap(),
|
|
895 |
/// The paramters can be specified using functions \ref lowerMap(),
|
808 |
896 |
/// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
|
809 |
|
/// \ref costMap(), \ref supplyMap() and \ref stSupply()
|
810 |
|
/// functions. For example,
|
|
897 |
/// \ref costMap(), \ref supplyMap(), \ref stSupply(),
|
|
898 |
/// \ref problemType(), \ref flowMap() and \ref potentialMap().
|
|
899 |
/// For example,
|
811 |
900 |
/// \code
|
812 |
901 |
/// NetworkSimplex<ListDigraph> ns(graph);
|
813 |
902 |
/// ns.boundMaps(lower, upper).costMap(cost)
|
814 |
903 |
/// .supplyMap(sup).run();
|
815 |
904 |
/// \endcode
|
816 |
905 |
///
|
817 |
906 |
/// This function can be called more than once. All the parameters
|
818 |
907 |
/// that have been given are kept for the next call, unless
|
819 |
908 |
/// \ref reset() is called, thus only the modified parameters
|
820 |
909 |
/// have to be set again. See \ref reset() for examples.
|
821 |
910 |
///
|
822 |
911 |
/// \param pivot_rule The pivot rule that will be used during the
|
823 |
912 |
/// algorithm. For more information see \ref PivotRule.
|
824 |
913 |
///
|
825 |
914 |
/// \return \c true if a feasible flow can be found.
|
826 |
915 |
bool run(PivotRule pivot_rule = BLOCK_SEARCH) {
|
827 |
916 |
return init() && start(pivot_rule);
|
828 |
917 |
}
|
829 |
918 |
|
830 |
919 |
/// \brief Reset all the parameters that have been given before.
|
831 |
920 |
///
|
832 |
921 |
/// This function resets all the paramaters that have been given
|
833 |
|
/// using \ref lowerMap(), \ref upperMap(), \ref capacityMap(),
|
834 |
|
/// \ref boundMaps(), \ref costMap(), \ref supplyMap() and
|
835 |
|
/// \ref stSupply() functions before.
|
|
922 |
/// before using functions \ref lowerMap(), \ref upperMap(),
|
|
923 |
/// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
|
|
924 |
/// \ref supplyMap(), \ref stSupply(), \ref problemType(),
|
|
925 |
/// \ref flowMap() and \ref potentialMap().
|
836 |
926 |
///
|
837 |
927 |
/// It is useful for multiple run() calls. If this function is not
|
838 |
928 |
/// used, all the parameters given before are kept for the next
|
839 |
929 |
/// \ref run() call.
|
840 |
930 |
///
|
841 |
931 |
/// For example,
|
842 |
932 |
/// \code
|
843 |
933 |
/// NetworkSimplex<ListDigraph> ns(graph);
|
844 |
934 |
///
|
845 |
935 |
/// // First run
|
846 |
936 |
/// ns.lowerMap(lower).capacityMap(cap).costMap(cost)
|
847 |
937 |
/// .supplyMap(sup).run();
|
848 |
938 |
///
|
849 |
939 |
/// // Run again with modified cost map (reset() is not called,
|
850 |
940 |
/// // so only the cost map have to be set again)
|
851 |
941 |
/// cost[e] += 100;
|
852 |
942 |
/// ns.costMap(cost).run();
|
853 |
943 |
///
|
854 |
944 |
/// // Run again from scratch using reset()
|
855 |
945 |
/// // (the lower bounds will be set to zero on all arcs)
|
856 |
946 |
/// ns.reset();
|
857 |
947 |
/// ns.capacityMap(cap).costMap(cost)
|
858 |
948 |
/// .supplyMap(sup).run();
|
859 |
949 |
/// \endcode
|
860 |
950 |
///
|
861 |
951 |
/// \return <tt>(*this)</tt>
|
862 |
952 |
NetworkSimplex& reset() {
|
863 |
953 |
delete _plower;
|
864 |
954 |
delete _pupper;
|
865 |
955 |
delete _pcost;
|
866 |
956 |
delete _psupply;
|
867 |
957 |
_plower = NULL;
|
868 |
958 |
_pupper = NULL;
|
869 |
959 |
_pcost = NULL;
|
870 |
960 |
_psupply = NULL;
|
871 |
961 |
_pstsup = false;
|
|
962 |
_ptype = GEQ;
|
|
963 |
if (_local_flow) delete _flow_map;
|
|
964 |
if (_local_potential) delete _potential_map;
|
|
965 |
_flow_map = NULL;
|
|
966 |
_potential_map = NULL;
|
|
967 |
_local_flow = false;
|
|
968 |
_local_potential = false;
|
|
969 |
|
872 |
970 |
return *this;
|
873 |
971 |
}
|
874 |
972 |
|
875 |
973 |
/// @}
|
876 |
974 |
|
877 |
975 |
/// \name Query Functions
|
878 |
976 |
/// The results of the algorithm can be obtained using these
|
879 |
977 |
/// functions.\n
|
880 |
978 |
/// The \ref run() function must be called before using them.
|
881 |
979 |
|
882 |
980 |
/// @{
|
883 |
981 |
|
884 |
982 |
/// \brief Return the total cost of the found flow.
|
885 |
983 |
///
|
886 |
984 |
/// This function returns the total cost of the found flow.
|
887 |
985 |
/// The complexity of the function is O(e).
|
888 |
986 |
///
|
889 |
987 |
/// \note The return type of the function can be specified as a
|
890 |
988 |
/// template parameter. For example,
|
891 |
989 |
/// \code
|
892 |
990 |
/// ns.totalCost<double>();
|
893 |
991 |
/// \endcode
|
894 |
992 |
/// It is useful if the total cost cannot be stored in the \c Cost
|
895 |
993 |
/// type of the algorithm, which is the default return type of the
|
896 |
994 |
/// function.
|
897 |
995 |
///
|
898 |
996 |
/// \pre \ref run() must be called before using this function.
|
899 |
997 |
template <typename Num>
|
900 |
998 |
Num totalCost() const {
|
901 |
999 |
Num c = 0;
|
902 |
1000 |
if (_pcost) {
|
903 |
1001 |
for (ArcIt e(_graph); e != INVALID; ++e)
|
904 |
1002 |
c += (*_flow_map)[e] * (*_pcost)[e];
|
905 |
1003 |
} else {
|
906 |
1004 |
for (ArcIt e(_graph); e != INVALID; ++e)
|
907 |
1005 |
c += (*_flow_map)[e];
|
908 |
1006 |
}
|
909 |
1007 |
return c;
|
910 |
1008 |
}
|
911 |
1009 |
|
912 |
1010 |
#ifndef DOXYGEN
|
913 |
1011 |
Cost totalCost() const {
|
914 |
1012 |
return totalCost<Cost>();
|
915 |
1013 |
}
|
916 |
1014 |
#endif
|
917 |
1015 |
|
918 |
1016 |
/// \brief Return the flow on the given arc.
|
919 |
1017 |
///
|
920 |
1018 |
/// This function returns the flow on the given arc.
|
921 |
1019 |
///
|
922 |
1020 |
/// \pre \ref run() must be called before using this function.
|
923 |
1021 |
Flow flow(const Arc& a) const {
|
924 |
1022 |
return (*_flow_map)[a];
|
925 |
1023 |
}
|
926 |
1024 |
|
927 |
1025 |
/// \brief Return a const reference to the flow map.
|
928 |
1026 |
///
|
929 |
1027 |
/// This function returns a const reference to an arc map storing
|
930 |
1028 |
/// the found flow.
|
931 |
1029 |
///
|
932 |
1030 |
/// \pre \ref run() must be called before using this function.
|
933 |
1031 |
const FlowMap& flowMap() const {
|
934 |
1032 |
return *_flow_map;
|
935 |
1033 |
}
|
936 |
1034 |
|
937 |
1035 |
/// \brief Return the potential (dual value) of the given node.
|
938 |
1036 |
///
|
939 |
1037 |
/// This function returns the potential (dual value) of the
|
940 |
1038 |
/// given node.
|
941 |
1039 |
///
|
942 |
1040 |
/// \pre \ref run() must be called before using this function.
|
943 |
1041 |
Cost potential(const Node& n) const {
|
944 |
1042 |
return (*_potential_map)[n];
|
945 |
1043 |
}
|
946 |
1044 |
|
947 |
1045 |
/// \brief Return a const reference to the potential map
|
948 |
1046 |
/// (the dual solution).
|
949 |
1047 |
///
|
950 |
1048 |
/// This function returns a const reference to a node map storing
|
951 |
1049 |
/// the found potentials, which form the dual solution of the
|
952 |
1050 |
/// \ref min_cost_flow "minimum cost flow" problem.
|
953 |
1051 |
///
|
954 |
1052 |
/// \pre \ref run() must be called before using this function.
|
955 |
1053 |
const PotentialMap& potentialMap() const {
|
956 |
1054 |
return *_potential_map;
|
957 |
1055 |
}
|
958 |
1056 |
|
959 |
1057 |
/// @}
|
960 |
1058 |
|
961 |
1059 |
private:
|
962 |
1060 |
|
963 |
1061 |
// Initialize internal data structures
|
964 |
1062 |
bool init() {
|
965 |
1063 |
// Initialize result maps
|
966 |
1064 |
if (!_flow_map) {
|
967 |
1065 |
_flow_map = new FlowMap(_graph);
|
968 |
1066 |
_local_flow = true;
|
969 |
1067 |
}
|
970 |
1068 |
if (!_potential_map) {
|
971 |
1069 |
_potential_map = new PotentialMap(_graph);
|
972 |
1070 |
_local_potential = true;
|
973 |
1071 |
}
|
974 |
1072 |
|
975 |
1073 |
// Initialize vectors
|
976 |
1074 |
_node_num = countNodes(_graph);
|
977 |
1075 |
_arc_num = countArcs(_graph);
|
978 |
1076 |
int all_node_num = _node_num + 1;
|
979 |
1077 |
int all_arc_num = _arc_num + _node_num;
|
980 |
1078 |
if (_node_num == 0) return false;
|
981 |
1079 |
|
982 |
1080 |
_arc_ref.resize(_arc_num);
|
983 |
1081 |
_source.resize(all_arc_num);
|
984 |
1082 |
_target.resize(all_arc_num);
|
985 |
1083 |
|
986 |
1084 |
_cap.resize(all_arc_num);
|
987 |
1085 |
_cost.resize(all_arc_num);
|
988 |
1086 |
_supply.resize(all_node_num);
|
989 |
1087 |
_flow.resize(all_arc_num);
|
990 |
1088 |
_pi.resize(all_node_num);
|
991 |
1089 |
|
992 |
1090 |
_parent.resize(all_node_num);
|
993 |
1091 |
_pred.resize(all_node_num);
|
994 |
1092 |
_forward.resize(all_node_num);
|
995 |
1093 |
_thread.resize(all_node_num);
|
996 |
1094 |
_rev_thread.resize(all_node_num);
|
997 |
1095 |
_succ_num.resize(all_node_num);
|
998 |
1096 |
_last_succ.resize(all_node_num);
|
999 |
1097 |
_state.resize(all_arc_num);
|
1000 |
1098 |
|
1001 |
1099 |
// Initialize node related data
|
1002 |
1100 |
bool valid_supply = true;
|
|
1101 |
Flow sum_supply = 0;
|
1003 |
1102 |
if (!_pstsup && !_psupply) {
|
1004 |
1103 |
_pstsup = true;
|
1005 |
1104 |
_psource = _ptarget = NodeIt(_graph);
|
1006 |
1105 |
_pstflow = 0;
|
1007 |
1106 |
}
|
1008 |
1107 |
if (_psupply) {
|
1009 |
|
Flow sum = 0;
|
1010 |
1108 |
int i = 0;
|
1011 |
1109 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
|
1012 |
1110 |
_node_id[n] = i;
|
1013 |
1111 |
_supply[i] = (*_psupply)[n];
|
1014 |
|
sum += _supply[i];
|
|
1112 |
sum_supply += _supply[i];
|
1015 |
1113 |
}
|
1016 |
|
valid_supply = (sum == 0);
|
|
1114 |
valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
|
|
1115 |
(_ptype == LEQ && sum_supply >= 0);
|
1017 |
1116 |
} else {
|
1018 |
1117 |
int i = 0;
|
1019 |
1118 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
|
1020 |
1119 |
_node_id[n] = i;
|
1021 |
1120 |
_supply[i] = 0;
|
1022 |
1121 |
}
|
1023 |
1122 |
_supply[_node_id[_psource]] = _pstflow;
|
1024 |
|
_supply[_node_id[_ptarget]] = -_pstflow;
|
|
1123 |
_supply[_node_id[_ptarget]] = -_pstflow;
|
1025 |
1124 |
}
|
1026 |
1125 |
if (!valid_supply) return false;
|
1027 |
1126 |
|
|
1127 |
// Infinite capacity value
|
|
1128 |
Flow inf_cap =
|
|
1129 |
std::numeric_limits<Flow>::has_infinity ?
|
|
1130 |
std::numeric_limits<Flow>::infinity() :
|
|
1131 |
std::numeric_limits<Flow>::max();
|
|
1132 |
|
|
1133 |
// Initialize artifical cost
|
|
1134 |
Cost art_cost;
|
|
1135 |
if (std::numeric_limits<Cost>::is_exact) {
|
|
1136 |
art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
|
|
1137 |
} else {
|
|
1138 |
art_cost = std::numeric_limits<Cost>::min();
|
|
1139 |
for (int i = 0; i != _arc_num; ++i) {
|
|
1140 |
if (_cost[i] > art_cost) art_cost = _cost[i];
|
|
1141 |
}
|
|
1142 |
art_cost = (art_cost + 1) * _node_num;
|
|
1143 |
}
|
|
1144 |
|
|
1145 |
// Run Circulation to check if a feasible solution exists
|
|
1146 |
typedef ConstMap<Arc, Flow> ConstArcMap;
|
|
1147 |
FlowNodeMap *csup = NULL;
|
|
1148 |
bool local_csup = false;
|
|
1149 |
if (_psupply) {
|
|
1150 |
csup = _psupply;
|
|
1151 |
} else {
|
|
1152 |
csup = new FlowNodeMap(_graph, 0);
|
|
1153 |
(*csup)[_psource] = _pstflow;
|
|
1154 |
(*csup)[_ptarget] = -_pstflow;
|
|
1155 |
local_csup = true;
|
|
1156 |
}
|
|
1157 |
bool circ_result = false;
|
|
1158 |
if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
|
|
1159 |
// GEQ problem type
|
|
1160 |
if (_plower) {
|
|
1161 |
if (_pupper) {
|
|
1162 |
Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
|
|
1163 |
circ(_graph, *_plower, *_pupper, *csup);
|
|
1164 |
circ_result = circ.run();
|
|
1165 |
} else {
|
|
1166 |
Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
|
|
1167 |
circ(_graph, *_plower, ConstArcMap(inf_cap), *csup);
|
|
1168 |
circ_result = circ.run();
|
|
1169 |
}
|
|
1170 |
} else {
|
|
1171 |
if (_pupper) {
|
|
1172 |
Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
|
|
1173 |
circ(_graph, ConstArcMap(0), *_pupper, *csup);
|
|
1174 |
circ_result = circ.run();
|
|
1175 |
} else {
|
|
1176 |
Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
|
|
1177 |
circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *csup);
|
|
1178 |
circ_result = circ.run();
|
|
1179 |
}
|
|
1180 |
}
|
|
1181 |
} else {
|
|
1182 |
// LEQ problem type
|
|
1183 |
typedef ReverseDigraph<const GR> RevGraph;
|
|
1184 |
typedef NegMap<FlowNodeMap> NegNodeMap;
|
|
1185 |
RevGraph rgraph(_graph);
|
|
1186 |
NegNodeMap neg_csup(*csup);
|
|
1187 |
if (_plower) {
|
|
1188 |
if (_pupper) {
|
|
1189 |
Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
|
|
1190 |
circ(rgraph, *_plower, *_pupper, neg_csup);
|
|
1191 |
circ_result = circ.run();
|
|
1192 |
} else {
|
|
1193 |
Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
|
|
1194 |
circ(rgraph, *_plower, ConstArcMap(inf_cap), neg_csup);
|
|
1195 |
circ_result = circ.run();
|
|
1196 |
}
|
|
1197 |
} else {
|
|
1198 |
if (_pupper) {
|
|
1199 |
Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
|
|
1200 |
circ(rgraph, ConstArcMap(0), *_pupper, neg_csup);
|
|
1201 |
circ_result = circ.run();
|
|
1202 |
} else {
|
|
1203 |
Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
|
|
1204 |
circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup);
|
|
1205 |
circ_result = circ.run();
|
|
1206 |
}
|
|
1207 |
}
|
|
1208 |
}
|
|
1209 |
if (local_csup) delete csup;
|
|
1210 |
if (!circ_result) return false;
|
|
1211 |
|
1028 |
1212 |
// Set data for the artificial root node
|
1029 |
1213 |
_root = _node_num;
|
1030 |
1214 |
_parent[_root] = -1;
|
1031 |
1215 |
_pred[_root] = -1;
|
1032 |
1216 |
_thread[_root] = 0;
|
1033 |
1217 |
_rev_thread[0] = _root;
|
1034 |
1218 |
_succ_num[_root] = all_node_num;
|
1035 |
1219 |
_last_succ[_root] = _root - 1;
|
1036 |
|
_supply[_root] = 0;
|
1037 |
|
_pi[_root] = 0;
|
|
1220 |
_supply[_root] = -sum_supply;
|
|
1221 |
if (sum_supply < 0) {
|
|
1222 |
_pi[_root] = -art_cost;
|
|
1223 |
} else {
|
|
1224 |
_pi[_root] = art_cost;
|
|
1225 |
}
|
1038 |
1226 |
|
1039 |
1227 |
// Store the arcs in a mixed order
|
1040 |
1228 |
int k = std::max(int(sqrt(_arc_num)), 10);
|
1041 |
1229 |
int i = 0;
|
1042 |
1230 |
for (ArcIt e(_graph); e != INVALID; ++e) {
|
1043 |
1231 |
_arc_ref[i] = e;
|
1044 |
1232 |
if ((i += k) >= _arc_num) i = (i % k) + 1;
|
1045 |
1233 |
}
|
1046 |
1234 |
|
1047 |
1235 |
// Initialize arc maps
|
1048 |
|
Flow inf_cap =
|
1049 |
|
std::numeric_limits<Flow>::has_infinity ?
|
1050 |
|
std::numeric_limits<Flow>::infinity() :
|
1051 |
|
std::numeric_limits<Flow>::max();
|
1052 |
1236 |
if (_pupper && _pcost) {
|
1053 |
1237 |
for (int i = 0; i != _arc_num; ++i) {
|
1054 |
1238 |
Arc e = _arc_ref[i];
|
1055 |
1239 |
_source[i] = _node_id[_graph.source(e)];
|
1056 |
1240 |
_target[i] = _node_id[_graph.target(e)];
|
1057 |
1241 |
_cap[i] = (*_pupper)[e];
|
1058 |
1242 |
_cost[i] = (*_pcost)[e];
|
1059 |
1243 |
_flow[i] = 0;
|
1060 |
1244 |
_state[i] = STATE_LOWER;
|
1061 |
1245 |
}
|
1062 |
1246 |
} else {
|
1063 |
1247 |
for (int i = 0; i != _arc_num; ++i) {
|
1064 |
1248 |
Arc e = _arc_ref[i];
|
1065 |
1249 |
_source[i] = _node_id[_graph.source(e)];
|
1066 |
1250 |
_target[i] = _node_id[_graph.target(e)];
|
1067 |
1251 |
_flow[i] = 0;
|
1068 |
1252 |
_state[i] = STATE_LOWER;
|
1069 |
1253 |
}
|
1070 |
1254 |
if (_pupper) {
|
1071 |
1255 |
for (int i = 0; i != _arc_num; ++i)
|
1072 |
1256 |
_cap[i] = (*_pupper)[_arc_ref[i]];
|
1073 |
1257 |
} else {
|
1074 |
1258 |
for (int i = 0; i != _arc_num; ++i)
|
1075 |
1259 |
_cap[i] = inf_cap;
|
1076 |
1260 |
}
|
1077 |
1261 |
if (_pcost) {
|
1078 |
1262 |
for (int i = 0; i != _arc_num; ++i)
|
1079 |
1263 |
_cost[i] = (*_pcost)[_arc_ref[i]];
|
1080 |
1264 |
} else {
|
1081 |
1265 |
for (int i = 0; i != _arc_num; ++i)
|
1082 |
1266 |
_cost[i] = 1;
|
1083 |
1267 |
}
|
1084 |
1268 |
}
|
1085 |
1269 |
|
1086 |
|
// Initialize artifical cost
|
1087 |
|
Cost art_cost;
|
1088 |
|
if (std::numeric_limits<Cost>::is_exact) {
|
1089 |
|
art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
|
1090 |
|
} else {
|
1091 |
|
art_cost = std::numeric_limits<Cost>::min();
|
1092 |
|
for (int i = 0; i != _arc_num; ++i) {
|
1093 |
|
if (_cost[i] > art_cost) art_cost = _cost[i];
|
1094 |
|
}
|
1095 |
|
art_cost = (art_cost + 1) * _node_num;
|
1096 |
|
}
|
1097 |
|
|
1098 |
1270 |
// Remove non-zero lower bounds
|
1099 |
1271 |
if (_plower) {
|
1100 |
1272 |
for (int i = 0; i != _arc_num; ++i) {
|
1101 |
1273 |
Flow c = (*_plower)[_arc_ref[i]];
|
1102 |
1274 |
if (c != 0) {
|
1103 |
1275 |
_cap[i] -= c;
|
1104 |
1276 |
_supply[_source[i]] -= c;
|
1105 |
1277 |
_supply[_target[i]] += c;
|
1106 |
1278 |
}
|
1107 |
1279 |
}
|
1108 |
1280 |
}
|
1109 |
1281 |
|
1110 |
1282 |
// Add artificial arcs and initialize the spanning tree data structure
|
1111 |
1283 |
for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
|
1112 |
1284 |
_thread[u] = u + 1;
|
1113 |
1285 |
_rev_thread[u + 1] = u;
|
1114 |
1286 |
_succ_num[u] = 1;
|
1115 |
1287 |
_last_succ[u] = u;
|
1116 |
1288 |
_parent[u] = _root;
|
1117 |
1289 |
_pred[u] = e;
|
1118 |
1290 |
_cost[e] = art_cost;
|
1119 |
1291 |
_cap[e] = inf_cap;
|
1120 |
1292 |
_state[e] = STATE_TREE;
|
1121 |
|
if (_supply[u] >= 0) {
|
|
1293 |
if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
|
1122 |
1294 |
_flow[e] = _supply[u];
|
1123 |
1295 |
_forward[u] = true;
|
1124 |
|
_pi[u] = -art_cost;
|
|
1296 |
_pi[u] = -art_cost + _pi[_root];
|
1125 |
1297 |
} else {
|
1126 |
1298 |
_flow[e] = -_supply[u];
|
1127 |
1299 |
_forward[u] = false;
|
1128 |
|
_pi[u] = art_cost;
|
|
1300 |
_pi[u] = art_cost + _pi[_root];
|
1129 |
1301 |
}
|
1130 |
1302 |
}
|
1131 |
1303 |
|
1132 |
1304 |
return true;
|
1133 |
1305 |
}
|
1134 |
1306 |
|
1135 |
1307 |
// Find the join node
|
1136 |
1308 |
void findJoinNode() {
|
1137 |
1309 |
int u = _source[in_arc];
|
1138 |
1310 |
int v = _target[in_arc];
|
1139 |
1311 |
while (u != v) {
|
1140 |
1312 |
if (_succ_num[u] < _succ_num[v]) {
|
1141 |
1313 |
u = _parent[u];
|
1142 |
1314 |
} else {
|
1143 |
1315 |
v = _parent[v];
|
1144 |
1316 |
}
|
1145 |
1317 |
}
|
1146 |
1318 |
join = u;
|
1147 |
1319 |
}
|
1148 |
1320 |
|
1149 |
1321 |
// Find the leaving arc of the cycle and returns true if the
|
1150 |
1322 |
// leaving arc is not the same as the entering arc
|
1151 |
1323 |
bool findLeavingArc() {
|
1152 |
1324 |
// Initialize first and second nodes according to the direction
|
1153 |
1325 |
// of the cycle
|
1154 |
1326 |
if (_state[in_arc] == STATE_LOWER) {
|
1155 |
1327 |
first = _source[in_arc];
|
1156 |
1328 |
second = _target[in_arc];
|
1157 |
1329 |
} else {
|
1158 |
1330 |
first = _target[in_arc];
|
1159 |
1331 |
second = _source[in_arc];
|
1160 |
1332 |
}
|
1161 |
1333 |
delta = _cap[in_arc];
|
1162 |
1334 |
int result = 0;
|
1163 |
1335 |
Flow d;
|
1164 |
1336 |
int e;
|
1165 |
1337 |
|
1166 |
1338 |
// Search the cycle along the path form the first node to the root
|
1167 |
1339 |
for (int u = first; u != join; u = _parent[u]) {
|
1168 |
1340 |
e = _pred[u];
|
1169 |
1341 |
d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
|
1170 |
1342 |
if (d < delta) {
|
1171 |
1343 |
delta = d;
|
1172 |
1344 |
u_out = u;
|
1173 |
1345 |
result = 1;
|
1174 |
1346 |
}
|
1175 |
1347 |
}
|
1176 |
1348 |
// Search the cycle along the path form the second node to the root
|
1177 |
1349 |
for (int u = second; u != join; u = _parent[u]) {
|
1178 |
1350 |
e = _pred[u];
|
1179 |
1351 |
d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
|
1180 |
1352 |
if (d <= delta) {
|
1181 |
1353 |
delta = d;
|
1182 |
1354 |
u_out = u;
|
1183 |
1355 |
result = 2;
|
1184 |
1356 |
}
|
1185 |
1357 |
}
|
1186 |
1358 |
|
1187 |
1359 |
if (result == 1) {
|
1188 |
1360 |
u_in = first;
|
1189 |
1361 |
v_in = second;
|
1190 |
1362 |
} else {
|
1191 |
1363 |
u_in = second;
|
1192 |
1364 |
v_in = first;
|
1193 |
1365 |
}
|
1194 |
1366 |
return result != 0;
|
1195 |
1367 |
}
|
1196 |
1368 |
|
1197 |
1369 |
// Change _flow and _state vectors
|
1198 |
1370 |
void changeFlow(bool change) {
|
1199 |
1371 |
// Augment along the cycle
|
1200 |
1372 |
if (delta > 0) {
|
1201 |
1373 |
Flow val = _state[in_arc] * delta;
|
1202 |
1374 |
_flow[in_arc] += val;
|
1203 |
1375 |
for (int u = _source[in_arc]; u != join; u = _parent[u]) {
|
1204 |
1376 |
_flow[_pred[u]] += _forward[u] ? -val : val;
|
1205 |
1377 |
}
|
1206 |
1378 |
for (int u = _target[in_arc]; u != join; u = _parent[u]) {
|
1207 |
1379 |
_flow[_pred[u]] += _forward[u] ? val : -val;
|
1208 |
1380 |
}
|
1209 |
1381 |
}
|
1210 |
1382 |
// Update the state of the entering and leaving arcs
|
1211 |
1383 |
if (change) {
|
1212 |
1384 |
_state[in_arc] = STATE_TREE;
|
1213 |
1385 |
_state[_pred[u_out]] =
|
1214 |
1386 |
(_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
|
1215 |
1387 |
} else {
|
1216 |
1388 |
_state[in_arc] = -_state[in_arc];
|
1217 |
1389 |
}
|
1218 |
1390 |
}
|
1219 |
1391 |
|
1220 |
1392 |
// Update the tree structure
|
1221 |
1393 |
void updateTreeStructure() {
|
1222 |
1394 |
int u, w;
|
1223 |
1395 |
int old_rev_thread = _rev_thread[u_out];
|
1224 |
1396 |
int old_succ_num = _succ_num[u_out];
|
... |
... |
@@ -1289,127 +1461,122 @@
|
1289 |
1461 |
while (u != u_in) {
|
1290 |
1462 |
w = _parent[u];
|
1291 |
1463 |
_pred[u] = _pred[w];
|
1292 |
1464 |
_forward[u] = !_forward[w];
|
1293 |
1465 |
tmp_sc += _succ_num[u] - _succ_num[w];
|
1294 |
1466 |
_succ_num[u] = tmp_sc;
|
1295 |
1467 |
_last_succ[w] = tmp_ls;
|
1296 |
1468 |
u = w;
|
1297 |
1469 |
}
|
1298 |
1470 |
_pred[u_in] = in_arc;
|
1299 |
1471 |
_forward[u_in] = (u_in == _source[in_arc]);
|
1300 |
1472 |
_succ_num[u_in] = old_succ_num;
|
1301 |
1473 |
|
1302 |
1474 |
// Set limits for updating _last_succ form v_in and v_out
|
1303 |
1475 |
// towards the root
|
1304 |
1476 |
int up_limit_in = -1;
|
1305 |
1477 |
int up_limit_out = -1;
|
1306 |
1478 |
if (_last_succ[join] == v_in) {
|
1307 |
1479 |
up_limit_out = join;
|
1308 |
1480 |
} else {
|
1309 |
1481 |
up_limit_in = join;
|
1310 |
1482 |
}
|
1311 |
1483 |
|
1312 |
1484 |
// Update _last_succ from v_in towards the root
|
1313 |
1485 |
for (u = v_in; u != up_limit_in && _last_succ[u] == v_in;
|
1314 |
1486 |
u = _parent[u]) {
|
1315 |
1487 |
_last_succ[u] = _last_succ[u_out];
|
1316 |
1488 |
}
|
1317 |
1489 |
// Update _last_succ from v_out towards the root
|
1318 |
1490 |
if (join != old_rev_thread && v_in != old_rev_thread) {
|
1319 |
1491 |
for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
|
1320 |
1492 |
u = _parent[u]) {
|
1321 |
1493 |
_last_succ[u] = old_rev_thread;
|
1322 |
1494 |
}
|
1323 |
1495 |
} else {
|
1324 |
1496 |
for (u = v_out; u != up_limit_out && _last_succ[u] == old_last_succ;
|
1325 |
1497 |
u = _parent[u]) {
|
1326 |
1498 |
_last_succ[u] = _last_succ[u_out];
|
1327 |
1499 |
}
|
1328 |
1500 |
}
|
1329 |
1501 |
|
1330 |
1502 |
// Update _succ_num from v_in to join
|
1331 |
1503 |
for (u = v_in; u != join; u = _parent[u]) {
|
1332 |
1504 |
_succ_num[u] += old_succ_num;
|
1333 |
1505 |
}
|
1334 |
1506 |
// Update _succ_num from v_out to join
|
1335 |
1507 |
for (u = v_out; u != join; u = _parent[u]) {
|
1336 |
1508 |
_succ_num[u] -= old_succ_num;
|
1337 |
1509 |
}
|
1338 |
1510 |
}
|
1339 |
1511 |
|
1340 |
1512 |
// Update potentials
|
1341 |
1513 |
void updatePotential() {
|
1342 |
1514 |
Cost sigma = _forward[u_in] ?
|
1343 |
1515 |
_pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
|
1344 |
1516 |
_pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
|
1345 |
1517 |
// Update potentials in the subtree, which has been moved
|
1346 |
1518 |
int end = _thread[_last_succ[u_in]];
|
1347 |
1519 |
for (int u = u_in; u != end; u = _thread[u]) {
|
1348 |
1520 |
_pi[u] += sigma;
|
1349 |
1521 |
}
|
1350 |
1522 |
}
|
1351 |
1523 |
|
1352 |
1524 |
// Execute the algorithm
|
1353 |
1525 |
bool start(PivotRule pivot_rule) {
|
1354 |
1526 |
// Select the pivot rule implementation
|
1355 |
1527 |
switch (pivot_rule) {
|
1356 |
1528 |
case FIRST_ELIGIBLE:
|
1357 |
1529 |
return start<FirstEligiblePivotRule>();
|
1358 |
1530 |
case BEST_ELIGIBLE:
|
1359 |
1531 |
return start<BestEligiblePivotRule>();
|
1360 |
1532 |
case BLOCK_SEARCH:
|
1361 |
1533 |
return start<BlockSearchPivotRule>();
|
1362 |
1534 |
case CANDIDATE_LIST:
|
1363 |
1535 |
return start<CandidateListPivotRule>();
|
1364 |
1536 |
case ALTERING_LIST:
|
1365 |
1537 |
return start<AlteringListPivotRule>();
|
1366 |
1538 |
}
|
1367 |
1539 |
return false;
|
1368 |
1540 |
}
|
1369 |
1541 |
|
1370 |
1542 |
template <typename PivotRuleImpl>
|
1371 |
1543 |
bool start() {
|
1372 |
1544 |
PivotRuleImpl pivot(*this);
|
1373 |
1545 |
|
1374 |
1546 |
// Execute the Network Simplex algorithm
|
1375 |
1547 |
while (pivot.findEnteringArc()) {
|
1376 |
1548 |
findJoinNode();
|
1377 |
1549 |
bool change = findLeavingArc();
|
1378 |
1550 |
changeFlow(change);
|
1379 |
1551 |
if (change) {
|
1380 |
1552 |
updateTreeStructure();
|
1381 |
1553 |
updatePotential();
|
1382 |
1554 |
}
|
1383 |
1555 |
}
|
1384 |
1556 |
|
1385 |
|
// Check if the flow amount equals zero on all the artificial arcs
|
1386 |
|
for (int e = _arc_num; e != _arc_num + _node_num; ++e) {
|
1387 |
|
if (_flow[e] > 0) return false;
|
1388 |
|
}
|
1389 |
|
|
1390 |
1557 |
// Copy flow values to _flow_map
|
1391 |
1558 |
if (_plower) {
|
1392 |
1559 |
for (int i = 0; i != _arc_num; ++i) {
|
1393 |
1560 |
Arc e = _arc_ref[i];
|
1394 |
1561 |
_flow_map->set(e, (*_plower)[e] + _flow[i]);
|
1395 |
1562 |
}
|
1396 |
1563 |
} else {
|
1397 |
1564 |
for (int i = 0; i != _arc_num; ++i) {
|
1398 |
1565 |
_flow_map->set(_arc_ref[i], _flow[i]);
|
1399 |
1566 |
}
|
1400 |
1567 |
}
|
1401 |
1568 |
// Copy potential values to _potential_map
|
1402 |
1569 |
for (NodeIt n(_graph); n != INVALID; ++n) {
|
1403 |
1570 |
_potential_map->set(n, _pi[_node_id[n]]);
|
1404 |
1571 |
}
|
1405 |
1572 |
|
1406 |
1573 |
return true;
|
1407 |
1574 |
}
|
1408 |
1575 |
|
1409 |
1576 |
}; //class NetworkSimplex
|
1410 |
1577 |
|
1411 |
1578 |
///@}
|
1412 |
1579 |
|
1413 |
1580 |
} //namespace lemon
|
1414 |
1581 |
|
1415 |
1582 |
#endif //LEMON_NETWORK_SIMPLEX_H
|