No children
0
2
0
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-2010 |
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_CAPACITY_SCALING_H |
20 | 20 |
#define LEMON_CAPACITY_SCALING_H |
21 | 21 |
|
22 | 22 |
/// \ingroup min_cost_flow_algs |
23 | 23 |
/// |
24 | 24 |
/// \file |
25 | 25 |
/// \brief Capacity Scaling algorithm for finding a minimum cost flow. |
26 | 26 |
|
27 | 27 |
#include <vector> |
28 | 28 |
#include <limits> |
29 | 29 |
#include <lemon/core.h> |
30 | 30 |
#include <lemon/bin_heap.h> |
31 | 31 |
|
32 | 32 |
namespace lemon { |
33 | 33 |
|
34 | 34 |
/// \brief Default traits class of CapacityScaling algorithm. |
35 | 35 |
/// |
36 | 36 |
/// Default traits class of CapacityScaling algorithm. |
37 | 37 |
/// \tparam GR Digraph type. |
38 | 38 |
/// \tparam V The number type used for flow amounts, capacity bounds |
39 | 39 |
/// and supply values. By default it is \c int. |
40 | 40 |
/// \tparam C The number type used for costs and potentials. |
41 | 41 |
/// By default it is the same as \c V. |
42 | 42 |
template <typename GR, typename V = int, typename C = V> |
43 | 43 |
struct CapacityScalingDefaultTraits |
44 | 44 |
{ |
45 | 45 |
/// The type of the digraph |
46 | 46 |
typedef GR Digraph; |
47 | 47 |
/// The type of the flow amounts, capacity bounds and supply values |
48 | 48 |
typedef V Value; |
49 | 49 |
/// The type of the arc costs |
50 | 50 |
typedef C Cost; |
51 | 51 |
|
52 | 52 |
/// \brief The type of the heap used for internal Dijkstra computations. |
53 | 53 |
/// |
54 | 54 |
/// The type of the heap used for internal Dijkstra computations. |
55 | 55 |
/// It must conform to the \ref lemon::concepts::Heap "Heap" concept, |
56 | 56 |
/// its priority type must be \c Cost and its cross reference type |
57 | 57 |
/// must be \ref RangeMap "RangeMap<int>". |
58 | 58 |
typedef BinHeap<Cost, RangeMap<int> > Heap; |
59 | 59 |
}; |
60 | 60 |
|
61 | 61 |
/// \addtogroup min_cost_flow_algs |
62 | 62 |
/// @{ |
63 | 63 |
|
64 | 64 |
/// \brief Implementation of the Capacity Scaling algorithm for |
65 | 65 |
/// finding a \ref min_cost_flow "minimum cost flow". |
66 | 66 |
/// |
67 | 67 |
/// \ref CapacityScaling implements the capacity scaling version |
68 | 68 |
/// of the successive shortest path algorithm for finding a |
69 | 69 |
/// \ref min_cost_flow "minimum cost flow" \ref amo93networkflows, |
70 | 70 |
/// \ref edmondskarp72theoretical. It is an efficient dual |
71 | 71 |
/// solution method. |
72 | 72 |
/// |
73 | 73 |
/// Most of the parameters of the problem (except for the digraph) |
74 | 74 |
/// can be given using separate functions, and the algorithm can be |
75 | 75 |
/// executed using the \ref run() function. If some parameters are not |
76 | 76 |
/// specified, then default values will be used. |
77 | 77 |
/// |
78 | 78 |
/// \tparam GR The digraph type the algorithm runs on. |
79 | 79 |
/// \tparam V The number type used for flow amounts, capacity bounds |
80 | 80 |
/// and supply values in the algorithm. By default, it is \c int. |
81 | 81 |
/// \tparam C The number type used for costs and potentials in the |
82 | 82 |
/// algorithm. By default, it is the same as \c V. |
83 | 83 |
/// \tparam TR The traits class that defines various types used by the |
84 | 84 |
/// algorithm. By default, it is \ref CapacityScalingDefaultTraits |
85 | 85 |
/// "CapacityScalingDefaultTraits<GR, V, C>". |
86 | 86 |
/// In most cases, this parameter should not be set directly, |
87 | 87 |
/// consider to use the named template parameters instead. |
88 | 88 |
/// |
89 | 89 |
/// \warning Both \c V and \c C must be signed number types. |
90 |
/// \warning All input data (capacities, supply values, and costs) must |
|
91 |
/// be integer. |
|
90 |
/// \warning Capacity bounds and supply values must be integer, but |
|
91 |
/// arc costs can be arbitrary real numbers. |
|
92 | 92 |
/// \warning This algorithm does not support negative costs for |
93 | 93 |
/// arcs having infinite upper bound. |
94 | 94 |
#ifdef DOXYGEN |
95 | 95 |
template <typename GR, typename V, typename C, typename TR> |
96 | 96 |
#else |
97 | 97 |
template < typename GR, typename V = int, typename C = V, |
98 | 98 |
typename TR = CapacityScalingDefaultTraits<GR, V, C> > |
99 | 99 |
#endif |
100 | 100 |
class CapacityScaling |
101 | 101 |
{ |
102 | 102 |
public: |
103 | 103 |
|
104 | 104 |
/// The type of the digraph |
105 | 105 |
typedef typename TR::Digraph Digraph; |
106 | 106 |
/// The type of the flow amounts, capacity bounds and supply values |
107 | 107 |
typedef typename TR::Value Value; |
108 | 108 |
/// The type of the arc costs |
109 | 109 |
typedef typename TR::Cost Cost; |
110 | 110 |
|
111 | 111 |
/// The type of the heap used for internal Dijkstra computations |
112 | 112 |
typedef typename TR::Heap Heap; |
113 | 113 |
|
114 | 114 |
/// The \ref CapacityScalingDefaultTraits "traits class" of the algorithm |
115 | 115 |
typedef TR Traits; |
116 | 116 |
|
117 | 117 |
public: |
118 | 118 |
|
119 | 119 |
/// \brief Problem type constants for the \c run() function. |
120 | 120 |
/// |
121 | 121 |
/// Enum type containing the problem type constants that can be |
122 | 122 |
/// returned by the \ref run() function of the algorithm. |
123 | 123 |
enum ProblemType { |
124 | 124 |
/// The problem has no feasible solution (flow). |
125 | 125 |
INFEASIBLE, |
126 | 126 |
/// The problem has optimal solution (i.e. it is feasible and |
127 | 127 |
/// bounded), and the algorithm has found optimal flow and node |
128 | 128 |
/// potentials (primal and dual solutions). |
129 | 129 |
OPTIMAL, |
130 | 130 |
/// The digraph contains an arc of negative cost and infinite |
131 | 131 |
/// upper bound. It means that the objective function is unbounded |
132 | 132 |
/// on that arc, however, note that it could actually be bounded |
133 | 133 |
/// over the feasible flows, but this algroithm cannot handle |
134 | 134 |
/// these cases. |
135 | 135 |
UNBOUNDED |
136 | 136 |
}; |
137 | 137 |
|
138 | 138 |
private: |
139 | 139 |
|
140 | 140 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
141 | 141 |
|
142 | 142 |
typedef std::vector<int> IntVector; |
143 | 143 |
typedef std::vector<Value> ValueVector; |
144 | 144 |
typedef std::vector<Cost> CostVector; |
145 | 145 |
typedef std::vector<char> BoolVector; |
146 | 146 |
// Note: vector<char> is used instead of vector<bool> for efficiency reasons |
147 | 147 |
|
148 | 148 |
private: |
149 | 149 |
|
150 | 150 |
// Data related to the underlying digraph |
151 | 151 |
const GR &_graph; |
152 | 152 |
int _node_num; |
153 | 153 |
int _arc_num; |
154 | 154 |
int _res_arc_num; |
155 | 155 |
int _root; |
156 | 156 |
|
157 | 157 |
// Parameters of the problem |
158 | 158 |
bool _have_lower; |
159 | 159 |
Value _sum_supply; |
160 | 160 |
|
161 | 161 |
// Data structures for storing the digraph |
162 | 162 |
IntNodeMap _node_id; |
163 | 163 |
IntArcMap _arc_idf; |
164 | 164 |
IntArcMap _arc_idb; |
165 | 165 |
IntVector _first_out; |
166 | 166 |
BoolVector _forward; |
167 | 167 |
IntVector _source; |
168 | 168 |
IntVector _target; |
169 | 169 |
IntVector _reverse; |
170 | 170 |
|
171 | 171 |
// Node and arc data |
172 | 172 |
ValueVector _lower; |
173 | 173 |
ValueVector _upper; |
174 | 174 |
CostVector _cost; |
175 | 175 |
ValueVector _supply; |
176 | 176 |
|
177 | 177 |
ValueVector _res_cap; |
178 | 178 |
CostVector _pi; |
179 | 179 |
ValueVector _excess; |
180 | 180 |
IntVector _excess_nodes; |
181 | 181 |
IntVector _deficit_nodes; |
182 | 182 |
|
183 | 183 |
Value _delta; |
184 | 184 |
int _factor; |
185 | 185 |
IntVector _pred; |
186 | 186 |
|
187 | 187 |
public: |
188 | 188 |
|
189 | 189 |
/// \brief Constant for infinite upper bounds (capacities). |
190 | 190 |
/// |
191 | 191 |
/// Constant for infinite upper bounds (capacities). |
192 | 192 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
193 | 193 |
/// \c std::numeric_limits<Value>::max() otherwise. |
194 | 194 |
const Value INF; |
195 | 195 |
|
196 | 196 |
private: |
197 | 197 |
|
198 | 198 |
// Special implementation of the Dijkstra algorithm for finding |
199 | 199 |
// shortest paths in the residual network of the digraph with |
200 | 200 |
// respect to the reduced arc costs and modifying the node |
201 | 201 |
// potentials according to the found distance labels. |
202 | 202 |
class ResidualDijkstra |
203 | 203 |
{ |
204 | 204 |
private: |
205 | 205 |
|
206 | 206 |
int _node_num; |
207 | 207 |
bool _geq; |
208 | 208 |
const IntVector &_first_out; |
209 | 209 |
const IntVector &_target; |
210 | 210 |
const CostVector &_cost; |
211 | 211 |
const ValueVector &_res_cap; |
212 | 212 |
const ValueVector &_excess; |
213 | 213 |
CostVector &_pi; |
214 | 214 |
IntVector &_pred; |
215 | 215 |
|
216 | 216 |
IntVector _proc_nodes; |
217 | 217 |
CostVector _dist; |
218 | 218 |
|
219 | 219 |
public: |
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-2010 |
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_algs |
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 | 33 |
|
34 | 34 |
namespace lemon { |
35 | 35 |
|
36 | 36 |
/// \addtogroup min_cost_flow_algs |
37 | 37 |
/// @{ |
38 | 38 |
|
39 | 39 |
/// \brief Implementation of the primal Network Simplex algorithm |
40 | 40 |
/// for finding a \ref min_cost_flow "minimum cost flow". |
41 | 41 |
/// |
42 | 42 |
/// \ref NetworkSimplex implements the primal Network Simplex algorithm |
43 | 43 |
/// for finding a \ref min_cost_flow "minimum cost flow" |
44 | 44 |
/// \ref amo93networkflows, \ref dantzig63linearprog, |
45 | 45 |
/// \ref kellyoneill91netsimplex. |
46 | 46 |
/// This algorithm is a highly efficient specialized version of the |
47 | 47 |
/// linear programming simplex method directly for the minimum cost |
48 | 48 |
/// flow problem. |
49 | 49 |
/// |
50 | 50 |
/// In general, \ref NetworkSimplex and \ref CostScaling are the fastest |
51 | 51 |
/// implementations available in LEMON for this problem. |
52 | 52 |
/// Furthermore, this class supports both directions of the supply/demand |
53 | 53 |
/// inequality constraints. For more information, see \ref SupplyType. |
54 | 54 |
/// |
55 | 55 |
/// Most of the parameters of the problem (except for the digraph) |
56 | 56 |
/// can be given using separate functions, and the algorithm can be |
57 | 57 |
/// executed using the \ref run() function. If some parameters are not |
58 | 58 |
/// specified, then default values will be used. |
59 | 59 |
/// |
60 | 60 |
/// \tparam GR The digraph type the algorithm runs on. |
61 | 61 |
/// \tparam V The number type used for flow amounts, capacity bounds |
62 | 62 |
/// and supply values in the algorithm. By default, it is \c int. |
63 | 63 |
/// \tparam C The number type used for costs and potentials in the |
64 | 64 |
/// algorithm. By default, it is the same as \c V. |
65 | 65 |
/// |
66 | 66 |
/// \warning Both \c V and \c C must be signed number types. |
67 | 67 |
/// \warning All input data (capacities, supply values, and costs) must |
68 | 68 |
/// be integer. |
69 | 69 |
/// |
70 | 70 |
/// \note %NetworkSimplex provides five different pivot rule |
71 | 71 |
/// implementations, from which the most efficient one is used |
72 | 72 |
/// by default. For more information, see \ref PivotRule. |
73 | 73 |
template <typename GR, typename V = int, typename C = V> |
74 | 74 |
class NetworkSimplex |
75 | 75 |
{ |
76 | 76 |
public: |
77 | 77 |
|
78 | 78 |
/// The type of the flow amounts, capacity bounds and supply values |
79 | 79 |
typedef V Value; |
80 | 80 |
/// The type of the arc costs |
81 | 81 |
typedef C Cost; |
82 | 82 |
|
83 | 83 |
public: |
84 | 84 |
|
85 | 85 |
/// \brief Problem type constants for the \c run() function. |
86 | 86 |
/// |
87 | 87 |
/// Enum type containing the problem type constants that can be |
88 | 88 |
/// returned by the \ref run() function of the algorithm. |
89 | 89 |
enum ProblemType { |
90 | 90 |
/// The problem has no feasible solution (flow). |
91 | 91 |
INFEASIBLE, |
92 | 92 |
/// The problem has optimal solution (i.e. it is feasible and |
93 | 93 |
/// bounded), and the algorithm has found optimal flow and node |
94 | 94 |
/// potentials (primal and dual solutions). |
95 | 95 |
OPTIMAL, |
96 | 96 |
/// The objective function of the problem is unbounded, i.e. |
97 | 97 |
/// there is a directed cycle having negative total cost and |
98 | 98 |
/// infinite upper bound. |
99 | 99 |
UNBOUNDED |
100 | 100 |
}; |
101 | 101 |
|
102 | 102 |
/// \brief Constants for selecting the type of the supply constraints. |
103 | 103 |
/// |
104 | 104 |
/// Enum type containing constants for selecting the supply type, |
105 | 105 |
/// i.e. the direction of the inequalities in the supply/demand |
106 | 106 |
/// constraints of the \ref min_cost_flow "minimum cost flow problem". |
107 | 107 |
/// |
108 | 108 |
/// The default supply type is \c GEQ, the \c LEQ type can be |
109 | 109 |
/// selected using \ref supplyType(). |
110 | 110 |
/// The equality form is a special case of both supply types. |
111 | 111 |
enum SupplyType { |
112 | 112 |
/// This option means that there are <em>"greater or equal"</em> |
113 | 113 |
/// supply/demand constraints in the definition of the problem. |
114 | 114 |
GEQ, |
115 | 115 |
/// This option means that there are <em>"less or equal"</em> |
116 | 116 |
/// supply/demand constraints in the definition of the problem. |
117 | 117 |
LEQ |
118 | 118 |
}; |
119 | 119 |
|
120 | 120 |
/// \brief Constants for selecting the pivot rule. |
121 | 121 |
/// |
122 | 122 |
/// Enum type containing constants for selecting the pivot rule for |
123 | 123 |
/// the \ref run() function. |
124 | 124 |
/// |
125 |
/// \ref NetworkSimplex provides five different pivot rule |
|
126 |
/// implementations that significantly affect the running time |
|
125 |
/// \ref NetworkSimplex provides five different implementations for |
|
126 |
/// the pivot strategy that significantly affects the running time |
|
127 | 127 |
/// of the algorithm. |
128 |
/// By default, \ref BLOCK_SEARCH "Block Search" is used, which |
|
129 |
/// turend out to be the most efficient and the most robust on various |
|
130 |
/// test inputs. |
|
131 |
/// However, another pivot rule can be selected using the \ref run() |
|
132 |
/// |
|
128 |
/// According to experimental tests conducted on various problem |
|
129 |
/// instances, \ref BLOCK_SEARCH "Block Search" and |
|
130 |
/// \ref ALTERING_LIST "Altering Candidate List" rules turned out |
|
131 |
/// to be the most efficient. |
|
132 |
/// Since \ref BLOCK_SEARCH "Block Search" is a simpler strategy that |
|
133 |
/// seemed to be slightly more robust, it is used by default. |
|
134 |
/// However, another pivot rule can easily be selected using the |
|
135 |
/// \ref run() function with the proper parameter. |
|
133 | 136 |
enum PivotRule { |
134 | 137 |
|
135 | 138 |
/// The \e First \e Eligible pivot rule. |
136 | 139 |
/// The next eligible arc is selected in a wraparound fashion |
137 | 140 |
/// in every iteration. |
138 | 141 |
FIRST_ELIGIBLE, |
139 | 142 |
|
140 | 143 |
/// The \e Best \e Eligible pivot rule. |
141 | 144 |
/// The best eligible arc is selected in every iteration. |
142 | 145 |
BEST_ELIGIBLE, |
143 | 146 |
|
144 | 147 |
/// The \e Block \e Search pivot rule. |
145 | 148 |
/// A specified number of arcs are examined in every iteration |
146 | 149 |
/// in a wraparound fashion and the best eligible arc is selected |
147 | 150 |
/// from this block. |
148 | 151 |
BLOCK_SEARCH, |
149 | 152 |
|
150 | 153 |
/// The \e Candidate \e List pivot rule. |
151 | 154 |
/// In a major iteration a candidate list is built from eligible arcs |
152 | 155 |
/// in a wraparound fashion and in the following minor iterations |
153 | 156 |
/// the best eligible arc is selected from this list. |
154 | 157 |
CANDIDATE_LIST, |
155 | 158 |
|
156 | 159 |
/// The \e Altering \e Candidate \e List pivot rule. |
157 | 160 |
/// It is a modified version of the Candidate List method. |
158 |
/// It keeps only |
|
161 |
/// It keeps only a few of the best eligible arcs from the former |
|
159 | 162 |
/// candidate list and extends this list in every iteration. |
160 | 163 |
ALTERING_LIST |
161 | 164 |
}; |
162 | 165 |
|
163 | 166 |
private: |
164 | 167 |
|
165 | 168 |
TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
166 | 169 |
|
167 | 170 |
typedef std::vector<int> IntVector; |
168 | 171 |
typedef std::vector<Value> ValueVector; |
169 | 172 |
typedef std::vector<Cost> CostVector; |
170 | 173 |
typedef std::vector<signed char> CharVector; |
171 | 174 |
// Note: vector<signed char> is used instead of vector<ArcState> and |
172 | 175 |
// vector<ArcDirection> for efficiency reasons |
173 | 176 |
|
174 | 177 |
// State constants for arcs |
175 | 178 |
enum ArcState { |
176 | 179 |
STATE_UPPER = -1, |
177 | 180 |
STATE_TREE = 0, |
178 | 181 |
STATE_LOWER = 1 |
179 | 182 |
}; |
180 | 183 |
|
181 | 184 |
// Direction constants for tree arcs |
182 | 185 |
enum ArcDirection { |
183 | 186 |
DIR_DOWN = -1, |
184 | 187 |
DIR_UP = 1 |
185 | 188 |
}; |
186 | 189 |
|
187 | 190 |
private: |
188 | 191 |
|
189 | 192 |
// Data related to the underlying digraph |
190 | 193 |
const GR &_graph; |
191 | 194 |
int _node_num; |
192 | 195 |
int _arc_num; |
193 | 196 |
int _all_arc_num; |
194 | 197 |
int _search_arc_num; |
195 | 198 |
|
196 | 199 |
// Parameters of the problem |
197 | 200 |
bool _have_lower; |
198 | 201 |
SupplyType _stype; |
199 | 202 |
Value _sum_supply; |
200 | 203 |
|
201 | 204 |
// Data structures for storing the digraph |
202 | 205 |
IntNodeMap _node_id; |
203 | 206 |
IntArcMap _arc_id; |
204 | 207 |
IntVector _source; |
205 | 208 |
IntVector _target; |
206 | 209 |
bool _arc_mixing; |
207 | 210 |
|
208 | 211 |
// Node and arc data |
209 | 212 |
ValueVector _lower; |
210 | 213 |
ValueVector _upper; |
211 | 214 |
ValueVector _cap; |
212 | 215 |
CostVector _cost; |
213 | 216 |
ValueVector _supply; |
214 | 217 |
ValueVector _flow; |
215 | 218 |
CostVector _pi; |
216 | 219 |
|
217 | 220 |
// Data for storing the spanning tree structure |
218 | 221 |
IntVector _parent; |
219 | 222 |
IntVector _pred; |
220 | 223 |
IntVector _thread; |
221 | 224 |
IntVector _rev_thread; |
222 | 225 |
IntVector _succ_num; |
223 | 226 |
IntVector _last_succ; |
224 | 227 |
CharVector _pred_dir; |
225 | 228 |
CharVector _state; |
226 | 229 |
IntVector _dirty_revs; |
227 | 230 |
int _root; |
228 | 231 |
|
229 | 232 |
// Temporary data used in the current pivot iteration |
230 | 233 |
int in_arc, join, u_in, v_in, u_out, v_out; |
231 | 234 |
Value delta; |
232 | 235 |
|
233 | 236 |
const Value MAX; |
234 | 237 |
|
235 | 238 |
public: |
236 | 239 |
|
237 | 240 |
/// \brief Constant for infinite upper bounds (capacities). |
238 | 241 |
/// |
239 | 242 |
/// Constant for infinite upper bounds (capacities). |
240 | 243 |
/// It is \c std::numeric_limits<Value>::infinity() if available, |
241 | 244 |
/// \c std::numeric_limits<Value>::max() otherwise. |
242 | 245 |
const Value INF; |
243 | 246 |
|
244 | 247 |
private: |
245 | 248 |
|
246 | 249 |
// Implementation of the First Eligible pivot rule |
247 | 250 |
class FirstEligiblePivotRule |
248 | 251 |
{ |
249 | 252 |
private: |
250 | 253 |
|
251 | 254 |
// References to the NetworkSimplex class |
252 | 255 |
const IntVector &_source; |
253 | 256 |
const IntVector &_target; |
254 | 257 |
const CostVector &_cost; |
255 | 258 |
const CharVector &_state; |
256 | 259 |
const CostVector &_pi; |
257 | 260 |
int &_in_arc; |
258 | 261 |
int _search_arc_num; |
259 | 262 |
|
260 | 263 |
// Pivot rule data |
261 | 264 |
int _next_arc; |
262 | 265 |
|
263 | 266 |
public: |
264 | 267 |
|
265 | 268 |
// Constructor |
266 | 269 |
FirstEligiblePivotRule(NetworkSimplex &ns) : |
267 | 270 |
_source(ns._source), _target(ns._target), |
268 | 271 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
269 | 272 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
270 | 273 |
_next_arc(0) |
271 | 274 |
{} |
272 | 275 |
|
273 | 276 |
// Find next entering arc |
274 | 277 |
bool findEnteringArc() { |
275 | 278 |
Cost c; |
276 | 279 |
for (int e = _next_arc; e != _search_arc_num; ++e) { |
277 | 280 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
278 | 281 |
if (c < 0) { |
279 | 282 |
_in_arc = e; |
280 | 283 |
_next_arc = e + 1; |
281 | 284 |
return true; |
282 | 285 |
} |
283 | 286 |
} |
284 | 287 |
for (int e = 0; e != _next_arc; ++e) { |
285 | 288 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
286 | 289 |
if (c < 0) { |
... | ... |
@@ -413,343 +416,343 @@ |
413 | 416 |
private: |
414 | 417 |
|
415 | 418 |
// References to the NetworkSimplex class |
416 | 419 |
const IntVector &_source; |
417 | 420 |
const IntVector &_target; |
418 | 421 |
const CostVector &_cost; |
419 | 422 |
const CharVector &_state; |
420 | 423 |
const CostVector &_pi; |
421 | 424 |
int &_in_arc; |
422 | 425 |
int _search_arc_num; |
423 | 426 |
|
424 | 427 |
// Pivot rule data |
425 | 428 |
IntVector _candidates; |
426 | 429 |
int _list_length, _minor_limit; |
427 | 430 |
int _curr_length, _minor_count; |
428 | 431 |
int _next_arc; |
429 | 432 |
|
430 | 433 |
public: |
431 | 434 |
|
432 | 435 |
/// Constructor |
433 | 436 |
CandidateListPivotRule(NetworkSimplex &ns) : |
434 | 437 |
_source(ns._source), _target(ns._target), |
435 | 438 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
436 | 439 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
437 | 440 |
_next_arc(0) |
438 | 441 |
{ |
439 | 442 |
// The main parameters of the pivot rule |
440 | 443 |
const double LIST_LENGTH_FACTOR = 0.25; |
441 | 444 |
const int MIN_LIST_LENGTH = 10; |
442 | 445 |
const double MINOR_LIMIT_FACTOR = 0.1; |
443 | 446 |
const int MIN_MINOR_LIMIT = 3; |
444 | 447 |
|
445 | 448 |
_list_length = std::max( int(LIST_LENGTH_FACTOR * |
446 | 449 |
std::sqrt(double(_search_arc_num))), |
447 | 450 |
MIN_LIST_LENGTH ); |
448 | 451 |
_minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length), |
449 | 452 |
MIN_MINOR_LIMIT ); |
450 | 453 |
_curr_length = _minor_count = 0; |
451 | 454 |
_candidates.resize(_list_length); |
452 | 455 |
} |
453 | 456 |
|
454 | 457 |
/// Find next entering arc |
455 | 458 |
bool findEnteringArc() { |
456 | 459 |
Cost min, c; |
457 | 460 |
int e; |
458 | 461 |
if (_curr_length > 0 && _minor_count < _minor_limit) { |
459 | 462 |
// Minor iteration: select the best eligible arc from the |
460 | 463 |
// current candidate list |
461 | 464 |
++_minor_count; |
462 | 465 |
min = 0; |
463 | 466 |
for (int i = 0; i < _curr_length; ++i) { |
464 | 467 |
e = _candidates[i]; |
465 | 468 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
466 | 469 |
if (c < min) { |
467 | 470 |
min = c; |
468 | 471 |
_in_arc = e; |
469 | 472 |
} |
470 | 473 |
else if (c >= 0) { |
471 | 474 |
_candidates[i--] = _candidates[--_curr_length]; |
472 | 475 |
} |
473 | 476 |
} |
474 | 477 |
if (min < 0) return true; |
475 | 478 |
} |
476 | 479 |
|
477 | 480 |
// Major iteration: build a new candidate list |
478 | 481 |
min = 0; |
479 | 482 |
_curr_length = 0; |
480 | 483 |
for (e = _next_arc; e != _search_arc_num; ++e) { |
481 | 484 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
482 | 485 |
if (c < 0) { |
483 | 486 |
_candidates[_curr_length++] = e; |
484 | 487 |
if (c < min) { |
485 | 488 |
min = c; |
486 | 489 |
_in_arc = e; |
487 | 490 |
} |
488 | 491 |
if (_curr_length == _list_length) goto search_end; |
489 | 492 |
} |
490 | 493 |
} |
491 | 494 |
for (e = 0; e != _next_arc; ++e) { |
492 | 495 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
493 | 496 |
if (c < 0) { |
494 | 497 |
_candidates[_curr_length++] = e; |
495 | 498 |
if (c < min) { |
496 | 499 |
min = c; |
497 | 500 |
_in_arc = e; |
498 | 501 |
} |
499 | 502 |
if (_curr_length == _list_length) goto search_end; |
500 | 503 |
} |
501 | 504 |
} |
502 | 505 |
if (_curr_length == 0) return false; |
503 | 506 |
|
504 | 507 |
search_end: |
505 | 508 |
_minor_count = 1; |
506 | 509 |
_next_arc = e; |
507 | 510 |
return true; |
508 | 511 |
} |
509 | 512 |
|
510 | 513 |
}; //class CandidateListPivotRule |
511 | 514 |
|
512 | 515 |
|
513 | 516 |
// Implementation of the Altering Candidate List pivot rule |
514 | 517 |
class AlteringListPivotRule |
515 | 518 |
{ |
516 | 519 |
private: |
517 | 520 |
|
518 | 521 |
// References to the NetworkSimplex class |
519 | 522 |
const IntVector &_source; |
520 | 523 |
const IntVector &_target; |
521 | 524 |
const CostVector &_cost; |
522 | 525 |
const CharVector &_state; |
523 | 526 |
const CostVector &_pi; |
524 | 527 |
int &_in_arc; |
525 | 528 |
int _search_arc_num; |
526 | 529 |
|
527 | 530 |
// Pivot rule data |
528 | 531 |
int _block_size, _head_length, _curr_length; |
529 | 532 |
int _next_arc; |
530 | 533 |
IntVector _candidates; |
531 | 534 |
CostVector _cand_cost; |
532 | 535 |
|
533 | 536 |
// Functor class to compare arcs during sort of the candidate list |
534 | 537 |
class SortFunc |
535 | 538 |
{ |
536 | 539 |
private: |
537 | 540 |
const CostVector &_map; |
538 | 541 |
public: |
539 | 542 |
SortFunc(const CostVector &map) : _map(map) {} |
540 | 543 |
bool operator()(int left, int right) { |
541 |
return _map[left] |
|
544 |
return _map[left] < _map[right]; |
|
542 | 545 |
} |
543 | 546 |
}; |
544 | 547 |
|
545 | 548 |
SortFunc _sort_func; |
546 | 549 |
|
547 | 550 |
public: |
548 | 551 |
|
549 | 552 |
// Constructor |
550 | 553 |
AlteringListPivotRule(NetworkSimplex &ns) : |
551 | 554 |
_source(ns._source), _target(ns._target), |
552 | 555 |
_cost(ns._cost), _state(ns._state), _pi(ns._pi), |
553 | 556 |
_in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num), |
554 | 557 |
_next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost) |
555 | 558 |
{ |
556 | 559 |
// The main parameters of the pivot rule |
557 | 560 |
const double BLOCK_SIZE_FACTOR = 1.0; |
558 | 561 |
const int MIN_BLOCK_SIZE = 10; |
559 |
const double HEAD_LENGTH_FACTOR = 0. |
|
562 |
const double HEAD_LENGTH_FACTOR = 0.01; |
|
560 | 563 |
const int MIN_HEAD_LENGTH = 3; |
561 | 564 |
|
562 | 565 |
_block_size = std::max( int(BLOCK_SIZE_FACTOR * |
563 | 566 |
std::sqrt(double(_search_arc_num))), |
564 | 567 |
MIN_BLOCK_SIZE ); |
565 | 568 |
_head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size), |
566 | 569 |
MIN_HEAD_LENGTH ); |
567 | 570 |
_candidates.resize(_head_length + _block_size); |
568 | 571 |
_curr_length = 0; |
569 | 572 |
} |
570 | 573 |
|
571 | 574 |
// Find next entering arc |
572 | 575 |
bool findEnteringArc() { |
573 | 576 |
// Check the current candidate list |
574 | 577 |
int e; |
575 | 578 |
Cost c; |
576 | 579 |
for (int i = 0; i != _curr_length; ++i) { |
577 | 580 |
e = _candidates[i]; |
578 | 581 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
579 | 582 |
if (c < 0) { |
580 | 583 |
_cand_cost[e] = c; |
581 | 584 |
} else { |
582 | 585 |
_candidates[i--] = _candidates[--_curr_length]; |
583 | 586 |
} |
584 | 587 |
} |
585 | 588 |
|
586 | 589 |
// Extend the list |
587 | 590 |
int cnt = _block_size; |
588 | 591 |
int limit = _head_length; |
589 | 592 |
|
590 | 593 |
for (e = _next_arc; e != _search_arc_num; ++e) { |
591 | 594 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
592 | 595 |
if (c < 0) { |
593 | 596 |
_cand_cost[e] = c; |
594 | 597 |
_candidates[_curr_length++] = e; |
595 | 598 |
} |
596 | 599 |
if (--cnt == 0) { |
597 | 600 |
if (_curr_length > limit) goto search_end; |
598 | 601 |
limit = 0; |
599 | 602 |
cnt = _block_size; |
600 | 603 |
} |
601 | 604 |
} |
602 | 605 |
for (e = 0; e != _next_arc; ++e) { |
603 |
_cand_cost[e] = _state[e] * |
|
604 |
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
|
605 |
|
|
606 |
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]); |
|
607 |
if (c < 0) { |
|
608 |
_cand_cost[e] = c; |
|
606 | 609 |
_candidates[_curr_length++] = e; |
607 | 610 |
} |
608 | 611 |
if (--cnt == 0) { |
609 | 612 |
if (_curr_length > limit) goto search_end; |
610 | 613 |
limit = 0; |
611 | 614 |
cnt = _block_size; |
612 | 615 |
} |
613 | 616 |
} |
614 | 617 |
if (_curr_length == 0) return false; |
615 | 618 |
|
616 | 619 |
search_end: |
617 | 620 |
|
618 |
// Make heap of the candidate list (approximating a partial sort) |
|
619 |
make_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
|
620 |
|
|
621 |
// Perform partial sort operation on the candidate list |
|
622 |
int new_length = std::min(_head_length + 1, _curr_length); |
|
623 |
std::partial_sort(_candidates.begin(), _candidates.begin() + new_length, |
|
624 |
_candidates.begin() + _curr_length, _sort_func); |
|
621 | 625 |
|
622 |
// |
|
626 |
// Select the entering arc and remove it from the list |
|
623 | 627 |
_in_arc = _candidates[0]; |
624 | 628 |
_next_arc = e; |
625 |
pop_heap( _candidates.begin(), _candidates.begin() + _curr_length, |
|
626 |
_sort_func ); |
|
627 |
|
|
629 |
_candidates[0] = _candidates[new_length - 1]; |
|
630 |
_curr_length = new_length - 1; |
|
628 | 631 |
return true; |
629 | 632 |
} |
630 | 633 |
|
631 | 634 |
}; //class AlteringListPivotRule |
632 | 635 |
|
633 | 636 |
public: |
634 | 637 |
|
635 | 638 |
/// \brief Constructor. |
636 | 639 |
/// |
637 | 640 |
/// The constructor of the class. |
638 | 641 |
/// |
639 | 642 |
/// \param graph The digraph the algorithm runs on. |
640 | 643 |
/// \param arc_mixing Indicate if the arcs will be stored in a |
641 | 644 |
/// mixed order in the internal data structure. |
642 | 645 |
/// In general, it leads to similar performance as using the original |
643 | 646 |
/// arc order, but it makes the algorithm more robust and in special |
644 | 647 |
/// cases, even significantly faster. Therefore, it is enabled by default. |
645 | 648 |
NetworkSimplex(const GR& graph, bool arc_mixing = true) : |
646 | 649 |
_graph(graph), _node_id(graph), _arc_id(graph), |
647 | 650 |
_arc_mixing(arc_mixing), |
648 | 651 |
MAX(std::numeric_limits<Value>::max()), |
649 | 652 |
INF(std::numeric_limits<Value>::has_infinity ? |
650 | 653 |
std::numeric_limits<Value>::infinity() : MAX) |
651 | 654 |
{ |
652 | 655 |
// Check the number types |
653 | 656 |
LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
654 | 657 |
"The flow type of NetworkSimplex must be signed"); |
655 | 658 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
656 | 659 |
"The cost type of NetworkSimplex must be signed"); |
657 | 660 |
|
658 | 661 |
// Reset data structures |
659 | 662 |
reset(); |
660 | 663 |
} |
661 | 664 |
|
662 | 665 |
/// \name Parameters |
663 | 666 |
/// The parameters of the algorithm can be specified using these |
664 | 667 |
/// functions. |
665 | 668 |
|
666 | 669 |
/// @{ |
667 | 670 |
|
668 | 671 |
/// \brief Set the lower bounds on the arcs. |
669 | 672 |
/// |
670 | 673 |
/// This function sets the lower bounds on the arcs. |
671 | 674 |
/// If it is not used before calling \ref run(), the lower bounds |
672 | 675 |
/// will be set to zero on all arcs. |
673 | 676 |
/// |
674 | 677 |
/// \param map An arc map storing the lower bounds. |
675 | 678 |
/// Its \c Value type must be convertible to the \c Value type |
676 | 679 |
/// of the algorithm. |
677 | 680 |
/// |
678 | 681 |
/// \return <tt>(*this)</tt> |
679 | 682 |
template <typename LowerMap> |
680 | 683 |
NetworkSimplex& lowerMap(const LowerMap& map) { |
681 | 684 |
_have_lower = true; |
682 | 685 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
683 | 686 |
_lower[_arc_id[a]] = map[a]; |
684 | 687 |
} |
685 | 688 |
return *this; |
686 | 689 |
} |
687 | 690 |
|
688 | 691 |
/// \brief Set the upper bounds (capacities) on the arcs. |
689 | 692 |
/// |
690 | 693 |
/// This function sets the upper bounds (capacities) on the arcs. |
691 | 694 |
/// If it is not used before calling \ref run(), the upper bounds |
692 | 695 |
/// will be set to \ref INF on all arcs (i.e. the flow value will be |
693 | 696 |
/// unbounded from above). |
694 | 697 |
/// |
695 | 698 |
/// \param map An arc map storing the upper bounds. |
696 | 699 |
/// Its \c Value type must be convertible to the \c Value type |
697 | 700 |
/// of the algorithm. |
698 | 701 |
/// |
699 | 702 |
/// \return <tt>(*this)</tt> |
700 | 703 |
template<typename UpperMap> |
701 | 704 |
NetworkSimplex& upperMap(const UpperMap& map) { |
702 | 705 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
703 | 706 |
_upper[_arc_id[a]] = map[a]; |
704 | 707 |
} |
705 | 708 |
return *this; |
706 | 709 |
} |
707 | 710 |
|
708 | 711 |
/// \brief Set the costs of the arcs. |
709 | 712 |
/// |
710 | 713 |
/// This function sets the costs of the arcs. |
711 | 714 |
/// If it is not used before calling \ref run(), the costs |
712 | 715 |
/// will be set to \c 1 on all arcs. |
713 | 716 |
/// |
714 | 717 |
/// \param map An arc map storing the costs. |
715 | 718 |
/// Its \c Value type must be convertible to the \c Cost type |
716 | 719 |
/// of the algorithm. |
717 | 720 |
/// |
718 | 721 |
/// \return <tt>(*this)</tt> |
719 | 722 |
template<typename CostMap> |
720 | 723 |
NetworkSimplex& costMap(const CostMap& map) { |
721 | 724 |
for (ArcIt a(_graph); a != INVALID; ++a) { |
722 | 725 |
_cost[_arc_id[a]] = map[a]; |
723 | 726 |
} |
724 | 727 |
return *this; |
725 | 728 |
} |
726 | 729 |
|
727 | 730 |
/// \brief Set the supply values of the nodes. |
728 | 731 |
/// |
729 | 732 |
/// This function sets the supply values of the nodes. |
730 | 733 |
/// If neither this function nor \ref stSupply() is used before |
731 | 734 |
/// calling \ref run(), the supply of each node will be set to zero. |
732 | 735 |
/// |
733 | 736 |
/// \param map A node map storing the supply values. |
734 | 737 |
/// Its \c Value type must be convertible to the \c Value type |
735 | 738 |
/// of the algorithm. |
736 | 739 |
/// |
737 | 740 |
/// \return <tt>(*this)</tt> |
738 | 741 |
/// |
739 | 742 |
/// \sa supplyType() |
740 | 743 |
template<typename SupplyMap> |
741 | 744 |
NetworkSimplex& supplyMap(const SupplyMap& map) { |
742 | 745 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
743 | 746 |
_supply[_node_id[n]] = map[n]; |
744 | 747 |
} |
745 | 748 |
return *this; |
746 | 749 |
} |
747 | 750 |
|
748 | 751 |
/// \brief Set single source and target nodes and a supply value. |
749 | 752 |
/// |
750 | 753 |
/// This function sets a single source node and a single target node |
751 | 754 |
/// and the required flow value. |
752 | 755 |
/// If neither this function nor \ref supplyMap() is used before |
753 | 756 |
/// calling \ref run(), the supply of each node will be set to zero. |
754 | 757 |
/// |
755 | 758 |
/// Using this function has the same effect as using \ref supplyMap() |
0 comments (0 inline)