1.1 --- a/lemon/network_simplex.h Sun Jan 13 10:26:55 2008 +0000
1.2 +++ b/lemon/network_simplex.h Sun Jan 13 10:32:14 2008 +0000
1.3 @@ -22,8 +22,7 @@
1.4 /// \ingroup min_cost_flow
1.5 ///
1.6 /// \file
1.7 -/// \brief The network simplex algorithm for finding a minimum cost
1.8 -/// flow.
1.9 +/// \brief The network simplex algorithm for finding a minimum cost flow.
1.10
1.11 #include <limits>
1.12 #include <lemon/graph_adaptor.h>
1.13 @@ -40,30 +39,29 @@
1.14 //#define _DEBUG_ITER_
1.15
1.16
1.17 -/// \brief State constant for edges at their lower bounds.
1.18 -#define LOWER 1
1.19 -/// \brief State constant for edges in the spanning tree.
1.20 -#define TREE 0
1.21 -/// \brief State constant for edges at their upper bounds.
1.22 -#define UPPER -1
1.23 +// State constant for edges at their lower bounds.
1.24 +#define LOWER 1
1.25 +// State constant for edges in the spanning tree.
1.26 +#define TREE 0
1.27 +// State constant for edges at their upper bounds.
1.28 +#define UPPER -1
1.29
1.30 #ifdef EDGE_BLOCK_PIVOT
1.31 #include <cmath>
1.32 - /// \brief Lower bound for the size of blocks.
1.33 - #define MIN_BLOCK_SIZE 10
1.34 + #define MIN_BLOCK_SIZE 10
1.35 #endif
1.36
1.37 #ifdef CANDIDATE_LIST_PIVOT
1.38 #include <vector>
1.39 - #define LIST_LENGTH_DIV 50
1.40 - #define MINOR_LIMIT_DIV 200
1.41 + #define LIST_LENGTH_DIV 50
1.42 + #define MINOR_LIMIT_DIV 200
1.43 #endif
1.44
1.45 #ifdef SORTED_LIST_PIVOT
1.46 #include <vector>
1.47 #include <algorithm>
1.48 #define LIST_LENGTH_DIV 100
1.49 - #define LOWER_DIV 4
1.50 + #define LOWER_DIV 4
1.51 #endif
1.52
1.53 namespace lemon {
1.54 @@ -74,8 +72,8 @@
1.55 /// \brief Implementation of the network simplex algorithm for
1.56 /// finding a minimum cost flow.
1.57 ///
1.58 - /// \ref lemon::NetworkSimplex "NetworkSimplex" implements the
1.59 - /// network simplex algorithm for finding a minimum cost flow.
1.60 + /// \ref NetworkSimplex implements the network simplex algorithm for
1.61 + /// finding a minimum cost flow.
1.62 ///
1.63 /// \param Graph The directed graph type the algorithm runs on.
1.64 /// \param LowerMap The type of the lower bound map.
1.65 @@ -84,12 +82,12 @@
1.66 /// \param SupplyMap The type of the supply map.
1.67 ///
1.68 /// \warning
1.69 - /// - Edge capacities and costs should be nonnegative integers.
1.70 - /// However \c CostMap::Value should be signed type.
1.71 + /// - Edge capacities and costs should be non-negative integers.
1.72 + /// However \c CostMap::Value should be signed type.
1.73 /// - Supply values should be signed integers.
1.74 /// - \c LowerMap::Value must be convertible to
1.75 - /// \c CapacityMap::Value and \c CapacityMap::Value must be
1.76 - /// convertible to \c SupplyMap::Value.
1.77 + /// \c CapacityMap::Value and \c CapacityMap::Value must be
1.78 + /// convertible to \c SupplyMap::Value.
1.79 ///
1.80 /// \author Peter Kovacs
1.81
1.82 @@ -107,12 +105,7 @@
1.83 typedef typename SupplyMap::Value Supply;
1.84
1.85 typedef SmartGraph SGraph;
1.86 - typedef typename SGraph::Node Node;
1.87 - typedef typename SGraph::NodeIt NodeIt;
1.88 - typedef typename SGraph::Edge Edge;
1.89 - typedef typename SGraph::EdgeIt EdgeIt;
1.90 - typedef typename SGraph::InEdgeIt InEdgeIt;
1.91 - typedef typename SGraph::OutEdgeIt OutEdgeIt;
1.92 + GRAPH_TYPEDEFS(typename SGraph);
1.93
1.94 typedef typename SGraph::template EdgeMap<Lower> SLowerMap;
1.95 typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
1.96 @@ -131,14 +124,14 @@
1.97
1.98 public:
1.99
1.100 - /// \brief The type of the flow map.
1.101 + /// The type of the flow map.
1.102 typedef typename Graph::template EdgeMap<Capacity> FlowMap;
1.103 - /// \brief The type of the potential map.
1.104 + /// The type of the potential map.
1.105 typedef typename Graph::template NodeMap<Cost> PotentialMap;
1.106
1.107 protected:
1.108
1.109 - /// \brief Map adaptor class for handling reduced edge costs.
1.110 + /// Map adaptor class for handling reduced edge costs.
1.111 class ReducedCostMap : public MapBase<Edge, Cost>
1.112 {
1.113 private:
1.114 @@ -150,65 +143,73 @@
1.115 public:
1.116
1.117 ReducedCostMap( const SGraph &_gr,
1.118 - const SCostMap &_cm,
1.119 - const SPotentialMap &_pm ) :
1.120 - gr(_gr), cost_map(_cm), pot_map(_pm) {}
1.121 + const SCostMap &_cm,
1.122 + const SPotentialMap &_pm ) :
1.123 + gr(_gr), cost_map(_cm), pot_map(_pm) {}
1.124
1.125 Cost operator[](const Edge &e) const {
1.126 - return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
1.127 + return cost_map[e] - pot_map[gr.source(e)] + pot_map[gr.target(e)];
1.128 }
1.129
1.130 }; //class ReducedCostMap
1.131
1.132 protected:
1.133
1.134 - /// \brief The directed graph the algorithm runs on.
1.135 + /// The directed graph the algorithm runs on.
1.136 SGraph graph;
1.137 - /// \brief The original graph.
1.138 + /// The original graph.
1.139 const Graph &graph_ref;
1.140 - /// \brief The original lower bound map.
1.141 + /// The original lower bound map.
1.142 const LowerMap *lower;
1.143 - /// \brief The capacity map.
1.144 + /// The capacity map.
1.145 SCapacityMap capacity;
1.146 - /// \brief The cost map.
1.147 + /// The cost map.
1.148 SCostMap cost;
1.149 - /// \brief The supply map.
1.150 + /// The supply map.
1.151 SSupplyMap supply;
1.152 - /// \brief The reduced cost map.
1.153 + /// The reduced cost map.
1.154 ReducedCostMap red_cost;
1.155 - /// \brief The sum of supply values equals zero.
1.156 bool valid_supply;
1.157
1.158 - /// \brief The edge map of the current flow.
1.159 + /// The edge map of the current flow.
1.160 SCapacityMap flow;
1.161 - /// \brief The edge map of the found flow on the original graph.
1.162 + /// The potential node map.
1.163 + SPotentialMap potential;
1.164 FlowMap flow_result;
1.165 - /// \brief The potential node map.
1.166 - SPotentialMap potential;
1.167 - /// \brief The potential node map on the original graph.
1.168 PotentialMap potential_result;
1.169
1.170 - /// \brief Node reference for the original graph.
1.171 + /// Node reference for the original graph.
1.172 NodeRefMap node_ref;
1.173 - /// \brief Edge reference for the original graph.
1.174 + /// Edge reference for the original graph.
1.175 EdgeRefMap edge_ref;
1.176
1.177 - /// \brief The depth node map of the spanning tree structure.
1.178 + /// The \c depth node map of the spanning tree structure.
1.179 IntNodeMap depth;
1.180 - /// \brief The parent node map of the spanning tree structure.
1.181 + /// The \c parent node map of the spanning tree structure.
1.182 NodeNodeMap parent;
1.183 - /// \brief The pred_edge node map of the spanning tree structure.
1.184 + /// The \c pred_edge node map of the spanning tree structure.
1.185 EdgeNodeMap pred_edge;
1.186 - /// \brief The thread node map of the spanning tree structure.
1.187 + /// The \c thread node map of the spanning tree structure.
1.188 NodeNodeMap thread;
1.189 - /// \brief The forward node map of the spanning tree structure.
1.190 + /// The \c forward node map of the spanning tree structure.
1.191 BoolNodeMap forward;
1.192 - /// \brief The state edge map.
1.193 + /// The state edge map.
1.194 IntEdgeMap state;
1.195 + /// The root node of the starting spanning tree.
1.196 + Node root;
1.197
1.198 + // The entering edge of the current pivot iteration.
1.199 + Edge in_edge;
1.200 + // Temporary nodes used in the current pivot iteration.
1.201 + Node join, u_in, v_in, u_out, v_out;
1.202 + Node right, first, second, last;
1.203 + Node stem, par_stem, new_stem;
1.204 + // The maximum augment amount along the found cycle in the current
1.205 + // pivot iteration.
1.206 + Capacity delta;
1.207
1.208 #ifdef EDGE_BLOCK_PIVOT
1.209 - /// \brief The size of blocks for the "Edge Block" pivot rule.
1.210 + /// The size of blocks for the "Edge Block" pivot rule.
1.211 int block_size;
1.212 #endif
1.213 #ifdef CANDIDATE_LIST_PIVOT
1.214 @@ -234,18 +235,6 @@
1.215 int list_index;
1.216 #endif
1.217
1.218 - // Root node of the starting spanning tree.
1.219 - Node root;
1.220 - // The entering edge of the current pivot iteration.
1.221 - Edge in_edge;
1.222 - // Temporary nodes used in the current pivot iteration.
1.223 - Node join, u_in, v_in, u_out, v_out;
1.224 - Node right, first, second, last;
1.225 - Node stem, par_stem, new_stem;
1.226 - // The maximum augment amount along the cycle in the current pivot
1.227 - // iteration.
1.228 - Capacity delta;
1.229 -
1.230 public :
1.231
1.232 /// \brief General constructor of the class (with lower bounds).
1.233 @@ -258,10 +247,10 @@
1.234 /// \param _cost The cost (length) values of the edges.
1.235 /// \param _supply The supply values of the nodes (signed).
1.236 NetworkSimplex( const Graph &_graph,
1.237 - const LowerMap &_lower,
1.238 - const CapacityMap &_capacity,
1.239 - const CostMap &_cost,
1.240 - const SupplyMap &_supply ) :
1.241 + const LowerMap &_lower,
1.242 + const CapacityMap &_capacity,
1.243 + const CostMap &_cost,
1.244 + const SupplyMap &_supply ) :
1.245 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
1.246 supply(graph), flow(graph), flow_result(_graph), potential(graph),
1.247 potential_result(_graph), depth(graph), parent(graph),
1.248 @@ -272,29 +261,29 @@
1.249 // Checking the sum of supply values
1.250 Supply sum = 0;
1.251 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
1.252 - sum += _supply[n];
1.253 + sum += _supply[n];
1.254 if (!(valid_supply = sum == 0)) return;
1.255
1.256 // Copying graph_ref to graph
1.257 graph.reserveNode(countNodes(graph_ref) + 1);
1.258 graph.reserveEdge(countEdges(graph_ref) + countNodes(graph_ref));
1.259 copyGraph(graph, graph_ref)
1.260 - .edgeMap(cost, _cost)
1.261 - .nodeRef(node_ref)
1.262 - .edgeRef(edge_ref)
1.263 - .run();
1.264 + .edgeMap(cost, _cost)
1.265 + .nodeRef(node_ref)
1.266 + .edgeRef(edge_ref)
1.267 + .run();
1.268
1.269 - // Removing nonzero lower bounds
1.270 + // Removing non-zero lower bounds
1.271 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
1.272 - capacity[edge_ref[e]] = _capacity[e] - _lower[e];
1.273 + capacity[edge_ref[e]] = _capacity[e] - _lower[e];
1.274 }
1.275 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
1.276 - Supply s = _supply[n];
1.277 - for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.278 - s += _lower[e];
1.279 - for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.280 - s -= _lower[e];
1.281 - supply[node_ref[n]] = s;
1.282 + Supply s = _supply[n];
1.283 + for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.284 + s += _lower[e];
1.285 + for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.286 + s -= _lower[e];
1.287 + supply[node_ref[n]] = s;
1.288 }
1.289 }
1.290
1.291 @@ -307,9 +296,9 @@
1.292 /// \param _cost The cost (length) values of the edges.
1.293 /// \param _supply The supply values of the nodes (signed).
1.294 NetworkSimplex( const Graph &_graph,
1.295 - const CapacityMap &_capacity,
1.296 - const CostMap &_cost,
1.297 - const SupplyMap &_supply ) :
1.298 + const CapacityMap &_capacity,
1.299 + const CostMap &_cost,
1.300 + const SupplyMap &_supply ) :
1.301 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
1.302 supply(graph), flow(graph), flow_result(_graph), potential(graph),
1.303 potential_result(_graph), depth(graph), parent(graph),
1.304 @@ -320,17 +309,17 @@
1.305 // Checking the sum of supply values
1.306 Supply sum = 0;
1.307 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
1.308 - sum += _supply[n];
1.309 + sum += _supply[n];
1.310 if (!(valid_supply = sum == 0)) return;
1.311
1.312 // Copying graph_ref to graph
1.313 copyGraph(graph, graph_ref)
1.314 - .edgeMap(capacity, _capacity)
1.315 - .edgeMap(cost, _cost)
1.316 - .nodeMap(supply, _supply)
1.317 - .nodeRef(node_ref)
1.318 - .edgeRef(edge_ref)
1.319 - .run();
1.320 + .edgeMap(capacity, _capacity)
1.321 + .edgeMap(cost, _cost)
1.322 + .nodeMap(supply, _supply)
1.323 + .nodeRef(node_ref)
1.324 + .edgeRef(edge_ref)
1.325 + .run();
1.326 }
1.327
1.328 /// \brief Simple constructor of the class (with lower bounds).
1.329 @@ -344,15 +333,14 @@
1.330 /// \param _s The source node.
1.331 /// \param _t The target node.
1.332 /// \param _flow_value The required amount of flow from node \c _s
1.333 - /// to node \c _t (i.e. the supply of \c _s and the demand of
1.334 - /// \c _t).
1.335 + /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
1.336 NetworkSimplex( const Graph &_graph,
1.337 - const LowerMap &_lower,
1.338 - const CapacityMap &_capacity,
1.339 - const CostMap &_cost,
1.340 - typename Graph::Node _s,
1.341 - typename Graph::Node _t,
1.342 - typename SupplyMap::Value _flow_value ) :
1.343 + const LowerMap &_lower,
1.344 + const CapacityMap &_capacity,
1.345 + const CostMap &_cost,
1.346 + typename Graph::Node _s,
1.347 + typename Graph::Node _t,
1.348 + typename SupplyMap::Value _flow_value ) :
1.349 graph_ref(_graph), lower(&_lower), capacity(graph), cost(graph),
1.350 supply(graph), flow(graph), flow_result(_graph), potential(graph),
1.351 potential_result(_graph), depth(graph), parent(graph),
1.352 @@ -362,24 +350,24 @@
1.353 {
1.354 // Copying graph_ref to graph
1.355 copyGraph(graph, graph_ref)
1.356 - .edgeMap(cost, _cost)
1.357 - .nodeRef(node_ref)
1.358 - .edgeRef(edge_ref)
1.359 - .run();
1.360 + .edgeMap(cost, _cost)
1.361 + .nodeRef(node_ref)
1.362 + .edgeRef(edge_ref)
1.363 + .run();
1.364
1.365 - // Removing nonzero lower bounds
1.366 + // Removing non-zero lower bounds
1.367 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e) {
1.368 - capacity[edge_ref[e]] = _capacity[e] - _lower[e];
1.369 + capacity[edge_ref[e]] = _capacity[e] - _lower[e];
1.370 }
1.371 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n) {
1.372 - Supply s = 0;
1.373 - if (n == _s) s = _flow_value;
1.374 - if (n == _t) s = -_flow_value;
1.375 - for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.376 - s += _lower[e];
1.377 - for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.378 - s -= _lower[e];
1.379 - supply[node_ref[n]] = s;
1.380 + Supply s = 0;
1.381 + if (n == _s) s = _flow_value;
1.382 + if (n == _t) s = -_flow_value;
1.383 + for (typename Graph::InEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.384 + s += _lower[e];
1.385 + for (typename Graph::OutEdgeIt e(graph_ref, n); e != INVALID; ++e)
1.386 + s -= _lower[e];
1.387 + supply[node_ref[n]] = s;
1.388 }
1.389 valid_supply = true;
1.390 }
1.391 @@ -394,14 +382,13 @@
1.392 /// \param _s The source node.
1.393 /// \param _t The target node.
1.394 /// \param _flow_value The required amount of flow from node \c _s
1.395 - /// to node \c _t (i.e. the supply of \c _s and the demand of
1.396 - /// \c _t).
1.397 + /// to node \c _t (i.e. the supply of \c _s and the demand of \c _t).
1.398 NetworkSimplex( const Graph &_graph,
1.399 - const CapacityMap &_capacity,
1.400 - const CostMap &_cost,
1.401 - typename Graph::Node _s,
1.402 - typename Graph::Node _t,
1.403 - typename SupplyMap::Value _flow_value ) :
1.404 + const CapacityMap &_capacity,
1.405 + const CostMap &_cost,
1.406 + typename Graph::Node _s,
1.407 + typename Graph::Node _t,
1.408 + typename SupplyMap::Value _flow_value ) :
1.409 graph_ref(_graph), lower(NULL), capacity(graph), cost(graph),
1.410 supply(graph, 0), flow(graph), flow_result(_graph), potential(graph),
1.411 potential_result(_graph), depth(graph), parent(graph),
1.412 @@ -411,16 +398,25 @@
1.413 {
1.414 // Copying graph_ref to graph
1.415 copyGraph(graph, graph_ref)
1.416 - .edgeMap(capacity, _capacity)
1.417 - .edgeMap(cost, _cost)
1.418 - .nodeRef(node_ref)
1.419 - .edgeRef(edge_ref)
1.420 - .run();
1.421 + .edgeMap(capacity, _capacity)
1.422 + .edgeMap(cost, _cost)
1.423 + .nodeRef(node_ref)
1.424 + .edgeRef(edge_ref)
1.425 + .run();
1.426 supply[node_ref[_s]] = _flow_value;
1.427 supply[node_ref[_t]] = -_flow_value;
1.428 valid_supply = true;
1.429 }
1.430
1.431 + /// \brief Runs the algorithm.
1.432 + ///
1.433 + /// Runs the algorithm.
1.434 + ///
1.435 + /// \return \c true if a feasible flow can be found.
1.436 + bool run() {
1.437 + return init() && start();
1.438 + }
1.439 +
1.440 /// \brief Returns a const reference to the flow map.
1.441 ///
1.442 /// Returns a const reference to the flow map.
1.443 @@ -450,19 +446,10 @@
1.444 Cost totalCost() const {
1.445 Cost c = 0;
1.446 for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
1.447 - c += flow_result[e] * cost[edge_ref[e]];
1.448 + c += flow_result[e] * cost[edge_ref[e]];
1.449 return c;
1.450 }
1.451
1.452 - /// \brief Runs the algorithm.
1.453 - ///
1.454 - /// Runs the algorithm.
1.455 - ///
1.456 - /// \return \c true if a feasible flow can be found.
1.457 - bool run() {
1.458 - return init() && start();
1.459 - }
1.460 -
1.461 protected:
1.462
1.463 /// \brief Extends the underlaying graph and initializes all the
1.464 @@ -472,8 +459,8 @@
1.465
1.466 // Initializing state and flow maps
1.467 for (EdgeIt e(graph); e != INVALID; ++e) {
1.468 - flow[e] = 0;
1.469 - state[e] = LOWER;
1.470 + flow[e] = 0;
1.471 + state[e] = LOWER;
1.472 }
1.473
1.474 // Adding an artificial root node to the graph
1.475 @@ -491,27 +478,27 @@
1.476 Edge e;
1.477 Cost max_cost = std::numeric_limits<Cost>::max() / 4;
1.478 for (NodeIt u(graph); u != INVALID; ++u) {
1.479 - if (u == root) continue;
1.480 - thread[last] = u;
1.481 - last = u;
1.482 - parent[u] = root;
1.483 - depth[u] = 1;
1.484 - sum += supply[u];
1.485 - if (supply[u] >= 0) {
1.486 - e = graph.addEdge(u, root);
1.487 - flow[e] = supply[u];
1.488 - forward[u] = true;
1.489 - potential[u] = max_cost;
1.490 - } else {
1.491 - e = graph.addEdge(root, u);
1.492 - flow[e] = -supply[u];
1.493 - forward[u] = false;
1.494 - potential[u] = -max_cost;
1.495 - }
1.496 - cost[e] = max_cost;
1.497 - capacity[e] = std::numeric_limits<Capacity>::max();
1.498 - state[e] = TREE;
1.499 - pred_edge[u] = e;
1.500 + if (u == root) continue;
1.501 + thread[last] = u;
1.502 + last = u;
1.503 + parent[u] = root;
1.504 + depth[u] = 1;
1.505 + sum += supply[u];
1.506 + if (supply[u] >= 0) {
1.507 + e = graph.addEdge(u, root);
1.508 + flow[e] = supply[u];
1.509 + forward[u] = true;
1.510 + potential[u] = max_cost;
1.511 + } else {
1.512 + e = graph.addEdge(root, u);
1.513 + flow[e] = -supply[u];
1.514 + forward[u] = false;
1.515 + potential[u] = -max_cost;
1.516 + }
1.517 + cost[e] = max_cost;
1.518 + capacity[e] = std::numeric_limits<Capacity>::max();
1.519 + state[e] = TREE;
1.520 + pred_edge[u] = e;
1.521 }
1.522 thread[last] = root;
1.523
1.524 @@ -520,8 +507,6 @@
1.525 int edge_num = countEdges(graph);
1.526 block_size = 2 * int(sqrt(countEdges(graph)));
1.527 if (block_size < MIN_BLOCK_SIZE) block_size = MIN_BLOCK_SIZE;
1.528 -// block_size = edge_num >= BLOCK_NUM * MIN_BLOCK_SIZE ?
1.529 -// edge_num / BLOCK_NUM : MIN_BLOCK_SIZE;
1.530 #endif
1.531 #ifdef CANDIDATE_LIST_PIVOT
1.532 int edge_num = countEdges(graph);
1.533 @@ -543,18 +528,18 @@
1.534 /// pivot rule.
1.535 bool findEnteringEdge(EdgeIt &next_edge) {
1.536 for (EdgeIt e = next_edge; e != INVALID; ++e) {
1.537 - if (state[e] * red_cost[e] < 0) {
1.538 - in_edge = e;
1.539 - next_edge = ++e;
1.540 - return true;
1.541 - }
1.542 + if (state[e] * red_cost[e] < 0) {
1.543 + in_edge = e;
1.544 + next_edge = ++e;
1.545 + return true;
1.546 + }
1.547 }
1.548 for (EdgeIt e(graph); e != next_edge; ++e) {
1.549 - if (state[e] * red_cost[e] < 0) {
1.550 - in_edge = e;
1.551 - next_edge = ++e;
1.552 - return true;
1.553 - }
1.554 + if (state[e] * red_cost[e] < 0) {
1.555 + in_edge = e;
1.556 + next_edge = ++e;
1.557 + return true;
1.558 + }
1.559 }
1.560 return false;
1.561 }
1.562 @@ -566,10 +551,10 @@
1.563 bool findEnteringEdge() {
1.564 Cost min = 0;
1.565 for (EdgeIt e(graph); e != INVALID; ++e) {
1.566 - if (state[e] * red_cost[e] < min) {
1.567 - min = state[e] * red_cost[e];
1.568 - in_edge = e;
1.569 - }
1.570 + if (state[e] * red_cost[e] < min) {
1.571 + min = state[e] * red_cost[e];
1.572 + in_edge = e;
1.573 + }
1.574 }
1.575 return min < 0;
1.576 }
1.577 @@ -584,30 +569,30 @@
1.578 EdgeIt min_edge(graph);
1.579 int cnt = 0;
1.580 for (EdgeIt e = next_edge; e != INVALID; ++e) {
1.581 - if ((curr = state[e] * red_cost[e]) < min) {
1.582 - min = curr;
1.583 - min_edge = e;
1.584 - }
1.585 - if (++cnt == block_size) {
1.586 - if (min < 0) break;
1.587 - cnt = 0;
1.588 - }
1.589 + if ((curr = state[e] * red_cost[e]) < min) {
1.590 + min = curr;
1.591 + min_edge = e;
1.592 + }
1.593 + if (++cnt == block_size) {
1.594 + if (min < 0) break;
1.595 + cnt = 0;
1.596 + }
1.597 }
1.598 if (!(min < 0)) {
1.599 - for (EdgeIt e(graph); e != next_edge; ++e) {
1.600 - if ((curr = state[e] * red_cost[e]) < min) {
1.601 - min = curr;
1.602 - min_edge = e;
1.603 - }
1.604 - if (++cnt == block_size) {
1.605 - if (min < 0) break;
1.606 - cnt = 0;
1.607 - }
1.608 - }
1.609 + for (EdgeIt e(graph); e != next_edge; ++e) {
1.610 + if ((curr = state[e] * red_cost[e]) < min) {
1.611 + min = curr;
1.612 + min_edge = e;
1.613 + }
1.614 + if (++cnt == block_size) {
1.615 + if (min < 0) break;
1.616 + cnt = 0;
1.617 + }
1.618 + }
1.619 }
1.620 in_edge = min_edge;
1.621 if ((next_edge = ++min_edge) == INVALID)
1.622 - next_edge = EdgeIt(graph);
1.623 + next_edge = EdgeIt(graph);
1.624 return min < 0;
1.625 }
1.626 #endif
1.627 @@ -619,15 +604,15 @@
1.628 typedef typename std::vector<Edge>::iterator ListIt;
1.629
1.630 if (minor_count >= minor_limit || candidates.size() == 0) {
1.631 - // Major iteration
1.632 - candidates.clear();
1.633 - for (EdgeIt e(graph); e != INVALID; ++e) {
1.634 - if (state[e] * red_cost[e] < 0) {
1.635 - candidates.push_back(e);
1.636 - if (candidates.size() == list_length) break;
1.637 - }
1.638 - }
1.639 - if (candidates.size() == 0) return false;
1.640 + // Major iteration
1.641 + candidates.clear();
1.642 + for (EdgeIt e(graph); e != INVALID; ++e) {
1.643 + if (state[e] * red_cost[e] < 0) {
1.644 + candidates.push_back(e);
1.645 + if (candidates.size() == list_length) break;
1.646 + }
1.647 + }
1.648 + if (candidates.size() == 0) return false;
1.649 }
1.650
1.651 // Minor iteration
1.652 @@ -636,10 +621,10 @@
1.653 Edge e;
1.654 for (int i = 0; i < candidates.size(); ++i) {
1.655 e = candidates[i];
1.656 - if (state[e] * red_cost[e] < min) {
1.657 - min = state[e] * red_cost[e];
1.658 - in_edge = e;
1.659 - }
1.660 + if (state[e] * red_cost[e] < min) {
1.661 + min = state[e] * red_cost[e];
1.662 + in_edge = e;
1.663 + }
1.664 }
1.665 return true;
1.666 }
1.667 @@ -655,9 +640,9 @@
1.668 const ReducedCostMap &rc;
1.669 public:
1.670 SortFunc(const IntEdgeMap &_st, const ReducedCostMap &_rc) :
1.671 - st(_st), rc(_rc) {}
1.672 + st(_st), rc(_rc) {}
1.673 bool operator()(const Edge &e1, const Edge &e2) {
1.674 - return st[e1] * rc[e1] < st[e2] * rc[e2];
1.675 + return st[e1] * rc[e1] < st[e2] * rc[e2];
1.676 }
1.677 };
1.678
1.679 @@ -668,19 +653,19 @@
1.680
1.681 // Minor iteration
1.682 while (list_index < candidates.size()) {
1.683 - in_edge = candidates[list_index++];
1.684 - if (state[in_edge] * red_cost[in_edge] < 0) return true;
1.685 + in_edge = candidates[list_index++];
1.686 + if (state[in_edge] * red_cost[in_edge] < 0) return true;
1.687 }
1.688
1.689 // Major iteration
1.690 candidates.clear();
1.691 Cost curr, min = 0;
1.692 for (EdgeIt e(graph); e != INVALID; ++e) {
1.693 - if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) {
1.694 - candidates.push_back(e);
1.695 - if (curr < min) min = curr;
1.696 - if (candidates.size() == list_length) break;
1.697 - }
1.698 + if ((curr = state[e] * red_cost[e]) < min / LOWER_DIV) {
1.699 + candidates.push_back(e);
1.700 + if (curr < min) min = curr;
1.701 + if (candidates.size() == list_length) break;
1.702 + }
1.703 }
1.704 if (candidates.size() == 0) return false;
1.705 sort(candidates.begin(), candidates.end(), sort_func);
1.706 @@ -696,12 +681,12 @@
1.707 Node u = graph.source(in_edge);
1.708 Node v = graph.target(in_edge);
1.709 while (u != v) {
1.710 - if (depth[u] == depth[v]) {
1.711 - u = parent[u];
1.712 - v = parent[v];
1.713 - }
1.714 - else if (depth[u] > depth[v]) u = parent[u];
1.715 - else v = parent[v];
1.716 + if (depth[u] == depth[v]) {
1.717 + u = parent[u];
1.718 + v = parent[v];
1.719 + }
1.720 + else if (depth[u] > depth[v]) u = parent[u];
1.721 + else v = parent[v];
1.722 }
1.723 return u;
1.724 }
1.725 @@ -712,11 +697,11 @@
1.726 // Initializing first and second nodes according to the direction
1.727 // of the cycle
1.728 if (state[in_edge] == LOWER) {
1.729 - first = graph.source(in_edge);
1.730 - second = graph.target(in_edge);
1.731 + first = graph.source(in_edge);
1.732 + second = graph.target(in_edge);
1.733 } else {
1.734 - first = graph.target(in_edge);
1.735 - second = graph.source(in_edge);
1.736 + first = graph.target(in_edge);
1.737 + second = graph.source(in_edge);
1.738 }
1.739 delta = capacity[in_edge];
1.740 bool result = false;
1.741 @@ -726,56 +711,56 @@
1.742 // Searching the cycle along the path form the first node to the
1.743 // root node
1.744 for (Node u = first; u != join; u = parent[u]) {
1.745 - e = pred_edge[u];
1.746 - d = forward[u] ? flow[e] : capacity[e] - flow[e];
1.747 - if (d < delta) {
1.748 - delta = d;
1.749 - u_out = u;
1.750 - u_in = first;
1.751 - v_in = second;
1.752 - result = true;
1.753 - }
1.754 + e = pred_edge[u];
1.755 + d = forward[u] ? flow[e] : capacity[e] - flow[e];
1.756 + if (d < delta) {
1.757 + delta = d;
1.758 + u_out = u;
1.759 + u_in = first;
1.760 + v_in = second;
1.761 + result = true;
1.762 + }
1.763 }
1.764 // Searching the cycle along the path form the second node to the
1.765 // root node
1.766 for (Node u = second; u != join; u = parent[u]) {
1.767 - e = pred_edge[u];
1.768 - d = forward[u] ? capacity[e] - flow[e] : flow[e];
1.769 - if (d <= delta) {
1.770 - delta = d;
1.771 - u_out = u;
1.772 - u_in = second;
1.773 - v_in = first;
1.774 - result = true;
1.775 - }
1.776 + e = pred_edge[u];
1.777 + d = forward[u] ? capacity[e] - flow[e] : flow[e];
1.778 + if (d <= delta) {
1.779 + delta = d;
1.780 + u_out = u;
1.781 + u_in = second;
1.782 + v_in = first;
1.783 + result = true;
1.784 + }
1.785 }
1.786 return result;
1.787 }
1.788
1.789 - /// \brief Changes flow and state edge maps.
1.790 + /// \brief Changes \c flow and \c state edge maps.
1.791 void changeFlows(bool change) {
1.792 // Augmenting along the cycle
1.793 if (delta > 0) {
1.794 - Capacity val = state[in_edge] * delta;
1.795 - flow[in_edge] += val;
1.796 - for (Node u = graph.source(in_edge); u != join; u = parent[u]) {
1.797 - flow[pred_edge[u]] += forward[u] ? -val : val;
1.798 - }
1.799 - for (Node u = graph.target(in_edge); u != join; u = parent[u]) {
1.800 - flow[pred_edge[u]] += forward[u] ? val : -val;
1.801 - }
1.802 + Capacity val = state[in_edge] * delta;
1.803 + flow[in_edge] += val;
1.804 + for (Node u = graph.source(in_edge); u != join; u = parent[u]) {
1.805 + flow[pred_edge[u]] += forward[u] ? -val : val;
1.806 + }
1.807 + for (Node u = graph.target(in_edge); u != join; u = parent[u]) {
1.808 + flow[pred_edge[u]] += forward[u] ? val : -val;
1.809 + }
1.810 }
1.811 // Updating the state of the entering and leaving edges
1.812 if (change) {
1.813 - state[in_edge] = TREE;
1.814 - state[pred_edge[u_out]] =
1.815 - (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER;
1.816 + state[in_edge] = TREE;
1.817 + state[pred_edge[u_out]] =
1.818 + (flow[pred_edge[u_out]] == 0) ? LOWER : UPPER;
1.819 } else {
1.820 - state[in_edge] = -state[in_edge];
1.821 + state[in_edge] = -state[in_edge];
1.822 }
1.823 }
1.824
1.825 - /// \brief Updates thread and parent node maps.
1.826 + /// \brief Updates \c thread and \c parent node maps.
1.827 void updateThreadParent() {
1.828 Node u;
1.829 v_out = parent[u_out];
1.830 @@ -783,12 +768,12 @@
1.831 // Handling the case when join and v_out coincide
1.832 bool par_first = false;
1.833 if (join == v_out) {
1.834 - for (u = join; u != u_in && u != v_in; u = thread[u]) ;
1.835 - if (u == v_in) {
1.836 - par_first = true;
1.837 - while (thread[u] != u_out) u = thread[u];
1.838 - first = u;
1.839 - }
1.840 + for (u = join; u != u_in && u != v_in; u = thread[u]) ;
1.841 + if (u == v_in) {
1.842 + par_first = true;
1.843 + while (thread[u] != u_out) u = thread[u];
1.844 + first = u;
1.845 + }
1.846 }
1.847
1.848 // Finding the last successor of u_in (u) and the node after it
1.849 @@ -796,8 +781,8 @@
1.850 for (u = u_in; depth[thread[u]] > depth[u_in]; u = thread[u]) ;
1.851 right = thread[u];
1.852 if (thread[v_in] == u_out) {
1.853 - for (last = u; depth[last] > depth[u_out]; last = thread[last]) ;
1.854 - if (last == u_out) last = thread[last];
1.855 + for (last = u; depth[last] > depth[u_out]; last = thread[last]) ;
1.856 + if (last == u_out) last = thread[last];
1.857 }
1.858 else last = thread[v_in];
1.859
1.860 @@ -805,65 +790,65 @@
1.861 thread[v_in] = stem = u_in;
1.862 par_stem = v_in;
1.863 while (stem != u_out) {
1.864 - thread[u] = new_stem = parent[stem];
1.865 + thread[u] = new_stem = parent[stem];
1.866
1.867 - // Finding the node just before the stem node (u) according to
1.868 - // the original thread index
1.869 - for (u = new_stem; thread[u] != stem; u = thread[u]) ;
1.870 - thread[u] = right;
1.871 + // Finding the node just before the stem node (u) according to
1.872 + // the original thread index
1.873 + for (u = new_stem; thread[u] != stem; u = thread[u]) ;
1.874 + thread[u] = right;
1.875
1.876 - // Changing the parent node of stem and shifting stem and
1.877 - // par_stem nodes
1.878 - parent[stem] = par_stem;
1.879 - par_stem = stem;
1.880 - stem = new_stem;
1.881 + // Changing the parent node of stem and shifting stem and
1.882 + // par_stem nodes
1.883 + parent[stem] = par_stem;
1.884 + par_stem = stem;
1.885 + stem = new_stem;
1.886
1.887 - // Finding the last successor of stem (u) and the node after it
1.888 - // (right) according to the thread index
1.889 - for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ;
1.890 - right = thread[u];
1.891 + // Finding the last successor of stem (u) and the node after it
1.892 + // (right) according to the thread index
1.893 + for (u = stem; depth[thread[u]] > depth[stem]; u = thread[u]) ;
1.894 + right = thread[u];
1.895 }
1.896 parent[u_out] = par_stem;
1.897 thread[u] = last;
1.898
1.899 if (join == v_out && par_first) {
1.900 - if (first != v_in) thread[first] = right;
1.901 + if (first != v_in) thread[first] = right;
1.902 } else {
1.903 - for (u = v_out; thread[u] != u_out; u = thread[u]) ;
1.904 - thread[u] = right;
1.905 + for (u = v_out; thread[u] != u_out; u = thread[u]) ;
1.906 + thread[u] = right;
1.907 }
1.908 }
1.909
1.910 - /// \brief Updates pred_edge and forward node maps.
1.911 + /// \brief Updates \c pred_edge and \c forward node maps.
1.912 void updatePredEdge() {
1.913 Node u = u_out, v;
1.914 while (u != u_in) {
1.915 - v = parent[u];
1.916 - pred_edge[u] = pred_edge[v];
1.917 - forward[u] = !forward[v];
1.918 - u = v;
1.919 + v = parent[u];
1.920 + pred_edge[u] = pred_edge[v];
1.921 + forward[u] = !forward[v];
1.922 + u = v;
1.923 }
1.924 pred_edge[u_in] = in_edge;
1.925 forward[u_in] = (u_in == graph.source(in_edge));
1.926 }
1.927
1.928 - /// \brief Updates depth and potential node maps.
1.929 + /// \brief Updates \c depth and \c potential node maps.
1.930 void updateDepthPotential() {
1.931 depth[u_in] = depth[v_in] + 1;
1.932 potential[u_in] = forward[u_in] ?
1.933 - potential[v_in] + cost[pred_edge[u_in]] :
1.934 - potential[v_in] - cost[pred_edge[u_in]];
1.935 + potential[v_in] + cost[pred_edge[u_in]] :
1.936 + potential[v_in] - cost[pred_edge[u_in]];
1.937
1.938 Node u = thread[u_in], v;
1.939 while (true) {
1.940 - v = parent[u];
1.941 - if (v == INVALID) break;
1.942 - depth[u] = depth[v] + 1;
1.943 - potential[u] = forward[u] ?
1.944 - potential[v] + cost[pred_edge[u]] :
1.945 - potential[v] - cost[pred_edge[u]];
1.946 - if (depth[u] <= depth[v_in]) break;
1.947 - u = thread[u];
1.948 + v = parent[u];
1.949 + if (v == INVALID) break;
1.950 + depth[u] = depth[v] + 1;
1.951 + potential[u] = forward[u] ?
1.952 + potential[v] + cost[pred_edge[u]] :
1.953 + potential[v] - cost[pred_edge[u]];
1.954 + if (depth[u] <= depth[v_in]) break;
1.955 + u = thread[u];
1.956 }
1.957 }
1.958
1.959 @@ -880,42 +865,42 @@
1.960 while (findEnteringEdge())
1.961 #endif
1.962 {
1.963 - join = findJoinNode();
1.964 - bool change = findLeavingEdge();
1.965 - changeFlows(change);
1.966 - if (change) {
1.967 - updateThreadParent();
1.968 - updatePredEdge();
1.969 - updateDepthPotential();
1.970 - }
1.971 + join = findJoinNode();
1.972 + bool change = findLeavingEdge();
1.973 + changeFlows(change);
1.974 + if (change) {
1.975 + updateThreadParent();
1.976 + updatePredEdge();
1.977 + updateDepthPotential();
1.978 + }
1.979 #ifdef _DEBUG_ITER_
1.980 - ++iter_num;
1.981 + ++iter_num;
1.982 #endif
1.983 }
1.984
1.985 #ifdef _DEBUG_ITER_
1.986 std::cout << "Network Simplex algorithm finished. " << iter_num
1.987 - << " pivot iterations performed." << std::endl;
1.988 + << " pivot iterations performed." << std::endl;
1.989 #endif
1.990
1.991 // Checking if the flow amount equals zero on all the
1.992 // artificial edges
1.993 for (InEdgeIt e(graph, root); e != INVALID; ++e)
1.994 - if (flow[e] > 0) return false;
1.995 + if (flow[e] > 0) return false;
1.996 for (OutEdgeIt e(graph, root); e != INVALID; ++e)
1.997 - if (flow[e] > 0) return false;
1.998 + if (flow[e] > 0) return false;
1.999
1.1000 // Copying flow values to flow_result
1.1001 if (lower) {
1.1002 - for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
1.1003 - flow_result[e] = (*lower)[e] + flow[edge_ref[e]];
1.1004 + for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
1.1005 + flow_result[e] = (*lower)[e] + flow[edge_ref[e]];
1.1006 } else {
1.1007 - for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
1.1008 - flow_result[e] = flow[edge_ref[e]];
1.1009 + for (typename Graph::EdgeIt e(graph_ref); e != INVALID; ++e)
1.1010 + flow_result[e] = flow[edge_ref[e]];
1.1011 }
1.1012 // Copying potential values to potential_result
1.1013 for (typename Graph::NodeIt n(graph_ref); n != INVALID; ++n)
1.1014 - potential_result[n] = potential[node_ref[n]];
1.1015 + potential_result[n] = potential[node_ref[n]];
1.1016
1.1017 return true;
1.1018 }