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