1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/lemon/cycle_canceling.h Mon May 07 11:42:18 2007 +0000
1.3 @@ -0,0 +1,509 @@
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_CYCLE_CANCELING_H
1.23 +#define LEMON_CYCLE_CANCELING_H
1.24 +
1.25 +/// \ingroup min_cost_flow
1.26 +///
1.27 +/// \file
1.28 +/// \brief A cycle canceling algorithm for finding a minimum cost flow.
1.29 +
1.30 +#include <vector>
1.31 +#include <lemon/circulation.h>
1.32 +#include <lemon/graph_adaptor.h>
1.33 +
1.34 +/// \brief The used cycle canceling method.
1.35 +#define LIMITED_CYCLE_CANCELING
1.36 +//#define MIN_MEAN_CYCLE_CANCELING
1.37 +
1.38 +#ifdef LIMITED_CYCLE_CANCELING
1.39 + #include <lemon/bellman_ford.h>
1.40 + /// \brief The maximum number of iterations for the first execution
1.41 + /// of the \ref lemon::BellmanFord "Bellman-Ford" algorithm.
1.42 + #define STARTING_LIMIT 2
1.43 + /// \brief The iteration limit for the
1.44 + /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by
1.45 + /// <tt>ALPHA_MUL % ALPHA_DIV</tt> in every round.
1.46 + #define ALPHA_MUL 3
1.47 + /// \brief The iteration limit for the
1.48 + /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by
1.49 + /// <tt>ALPHA_MUL % ALPHA_DIV</tt> in every round.
1.50 + #define ALPHA_DIV 2
1.51 +
1.52 +//#define _ONLY_ONE_CYCLE_
1.53 +//#define _DEBUG_ITER_
1.54 +#endif
1.55 +
1.56 +#ifdef MIN_MEAN_CYCLE_CANCELING
1.57 + #include <lemon/min_mean_cycle.h>
1.58 + #include <lemon/path.h>
1.59 +#endif
1.60 +
1.61 +namespace lemon {
1.62 +
1.63 + /// \addtogroup min_cost_flow
1.64 + /// @{
1.65 +
1.66 + /// \brief Implementation of a cycle canceling algorithm for finding
1.67 + /// a minimum cost flow.
1.68 + ///
1.69 + /// \ref lemon::CycleCanceling "CycleCanceling" implements a cycle
1.70 + /// canceling algorithm for finding a minimum cost flow.
1.71 + ///
1.72 + /// \param Graph The directed graph type the algorithm runs on.
1.73 + /// \param LowerMap The type of the lower bound map.
1.74 + /// \param CapacityMap The type of the capacity (upper bound) map.
1.75 + /// \param CostMap The type of the cost (length) map.
1.76 + /// \param SupplyMap The type of the supply map.
1.77 + ///
1.78 + /// \warning
1.79 + /// - Edge capacities and costs should be nonnegative integers.
1.80 + /// However \c CostMap::Value should be signed type.
1.81 + /// - Supply values should be integers.
1.82 + /// - \c LowerMap::Value must be convertible to
1.83 + /// \c CapacityMap::Value and \c CapacityMap::Value must be
1.84 + /// convertible to \c SupplyMap::Value.
1.85 + ///
1.86 + /// \author Peter Kovacs
1.87 +
1.88 +template < typename Graph,
1.89 + typename LowerMap = typename Graph::template EdgeMap<int>,
1.90 + typename CapacityMap = LowerMap,
1.91 + typename CostMap = typename Graph::template EdgeMap<int>,
1.92 + typename SupplyMap = typename Graph::template NodeMap
1.93 + <typename CapacityMap::Value> >
1.94 + class CycleCanceling
1.95 + {
1.96 + typedef typename Graph::Node Node;
1.97 + typedef typename Graph::NodeIt NodeIt;
1.98 + typedef typename Graph::Edge Edge;
1.99 + typedef typename Graph::EdgeIt EdgeIt;
1.100 + typedef typename Graph::InEdgeIt InEdgeIt;
1.101 + typedef typename Graph::OutEdgeIt OutEdgeIt;
1.102 +
1.103 + typedef typename LowerMap::Value Lower;
1.104 + typedef typename CapacityMap::Value Capacity;
1.105 + typedef typename CostMap::Value Cost;
1.106 + typedef typename SupplyMap::Value Supply;
1.107 + typedef typename Graph::template EdgeMap<Capacity> CapacityRefMap;
1.108 + typedef typename Graph::template NodeMap<Supply> SupplyRefMap;
1.109 +
1.110 + typedef ResGraphAdaptor< const Graph, Capacity,
1.111 + CapacityRefMap, CapacityRefMap > ResGraph;
1.112 + typedef typename ResGraph::Node ResNode;
1.113 + typedef typename ResGraph::NodeIt ResNodeIt;
1.114 + typedef typename ResGraph::Edge ResEdge;
1.115 + typedef typename ResGraph::EdgeIt ResEdgeIt;
1.116 +
1.117 + public:
1.118 +
1.119 + /// \brief The type of the flow map.
1.120 + typedef CapacityRefMap FlowMap;
1.121 +
1.122 + protected:
1.123 +
1.124 + /// \brief Map adaptor class for demand map.
1.125 + class DemandMap : public MapBase<Node, Supply>
1.126 + {
1.127 + private:
1.128 +
1.129 + const SupplyMap *map;
1.130 +
1.131 + public:
1.132 +
1.133 + typedef typename MapBase<Node, Supply>::Value Value;
1.134 + typedef typename MapBase<Node, Supply>::Key Key;
1.135 +
1.136 + DemandMap(const SupplyMap &_map) : map(&_map) {}
1.137 +
1.138 + Value operator[](const Key &e) const {
1.139 + return -(*map)[e];
1.140 + }
1.141 +
1.142 + }; //class DemandMap
1.143 +
1.144 + /// \brief Map adaptor class for handling residual edge costs.
1.145 + class ResCostMap : public MapBase<ResEdge, Cost>
1.146 + {
1.147 + private:
1.148 +
1.149 + const CostMap &cost_map;
1.150 +
1.151 + public:
1.152 +
1.153 + typedef typename MapBase<ResEdge, Cost>::Value Value;
1.154 + typedef typename MapBase<ResEdge, Cost>::Key Key;
1.155 +
1.156 + ResCostMap(const CostMap &_cost) : cost_map(_cost) {}
1.157 +
1.158 + Value operator[](const Key &e) const {
1.159 + return ResGraph::forward(e) ? cost_map[e] : -cost_map[e];
1.160 + }
1.161 +
1.162 + }; //class ResCostMap
1.163 +
1.164 + protected:
1.165 +
1.166 + /// \brief The directed graph the algorithm runs on.
1.167 + const Graph &graph;
1.168 + /// \brief The original lower bound map.
1.169 + const LowerMap *lower;
1.170 + /// \brief The modified capacity map.
1.171 + CapacityRefMap capacity;
1.172 + /// \brief The cost map.
1.173 + const CostMap &cost;
1.174 + /// \brief The modified supply map.
1.175 + SupplyRefMap supply;
1.176 + /// \brief The modified demand map.
1.177 + DemandMap demand;
1.178 + /// \brief The sum of supply values equals zero.
1.179 + bool valid_supply;
1.180 +
1.181 + /// \brief The current flow.
1.182 + FlowMap flow;
1.183 + /// \brief The residual graph.
1.184 + ResGraph res_graph;
1.185 + /// \brief The residual cost map.
1.186 + ResCostMap res_cost;
1.187 +
1.188 + public :
1.189 +
1.190 + /// \brief General constructor of the class (with lower bounds).
1.191 + ///
1.192 + /// General constructor of the class (with lower bounds).
1.193 + ///
1.194 + /// \param _graph The directed graph the algorithm runs on.
1.195 + /// \param _lower The lower bounds of the edges.
1.196 + /// \param _capacity The capacities (upper bounds) of the edges.
1.197 + /// \param _cost The cost (length) values of the edges.
1.198 + /// \param _supply The supply values of the nodes (signed).
1.199 + CycleCanceling( const Graph &_graph,
1.200 + const LowerMap &_lower,
1.201 + const CapacityMap &_capacity,
1.202 + const CostMap &_cost,
1.203 + const SupplyMap &_supply ) :
1.204 + graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
1.205 + supply(_graph), demand(supply), flow(_graph, 0),
1.206 + res_graph(_graph, capacity, flow), res_cost(cost)
1.207 + {
1.208 + // Removing nonzero lower bounds
1.209 + capacity = subMap(_capacity, _lower);
1.210 + Supply sum = 0;
1.211 + for (NodeIt n(graph); n != INVALID; ++n) {
1.212 + Supply s = _supply[n];
1.213 + for (InEdgeIt e(graph, n); e != INVALID; ++e)
1.214 + s += _lower[e];
1.215 + for (OutEdgeIt e(graph, n); e != INVALID; ++e)
1.216 + s -= _lower[e];
1.217 + sum += (supply[n] = s);
1.218 + }
1.219 + valid_supply = sum == 0;
1.220 + }
1.221 +
1.222 + /// \brief General constructor of the class (without lower bounds).
1.223 + ///
1.224 + /// General constructor of the class (without lower bounds).
1.225 + ///
1.226 + /// \param _graph The directed graph the algorithm runs on.
1.227 + /// \param _capacity The capacities (upper bounds) of the edges.
1.228 + /// \param _cost The cost (length) values of the edges.
1.229 + /// \param _supply The supply values of the nodes (signed).
1.230 + CycleCanceling( const Graph &_graph,
1.231 + const CapacityMap &_capacity,
1.232 + const CostMap &_cost,
1.233 + const SupplyMap &_supply ) :
1.234 + graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
1.235 + supply(_supply), demand(supply), flow(_graph, 0),
1.236 + res_graph(_graph, capacity, flow), res_cost(cost)
1.237 + {
1.238 + // Checking the sum of supply values
1.239 + Supply sum = 0;
1.240 + for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
1.241 + valid_supply = sum == 0;
1.242 + }
1.243 +
1.244 +
1.245 + /// \brief Simple constructor of the class (with lower bounds).
1.246 + ///
1.247 + /// Simple constructor of the class (with lower bounds).
1.248 + ///
1.249 + /// \param _graph The directed graph the algorithm runs on.
1.250 + /// \param _lower The lower bounds of the edges.
1.251 + /// \param _capacity The capacities (upper bounds) of the edges.
1.252 + /// \param _cost The cost (length) values of the edges.
1.253 + /// \param _s The source node.
1.254 + /// \param _t The target node.
1.255 + /// \param _flow_value The required amount of flow from node \c _s
1.256 + /// to node \c _t (i.e. the supply of \c _s and the demand of
1.257 + /// \c _t).
1.258 + CycleCanceling( const Graph &_graph,
1.259 + const LowerMap &_lower,
1.260 + const CapacityMap &_capacity,
1.261 + const CostMap &_cost,
1.262 + Node _s, Node _t,
1.263 + Supply _flow_value ) :
1.264 + graph(_graph), lower(&_lower), capacity(_graph), cost(_cost),
1.265 + supply(_graph), demand(supply), flow(_graph, 0),
1.266 + res_graph(_graph, capacity, flow), res_cost(cost)
1.267 + {
1.268 + // Removing nonzero lower bounds
1.269 + capacity = subMap(_capacity, _lower);
1.270 + for (NodeIt n(graph); n != INVALID; ++n) {
1.271 + Supply s = 0;
1.272 + if (n == _s) s = _flow_value;
1.273 + if (n == _t) s = -_flow_value;
1.274 + for (InEdgeIt e(graph, n); e != INVALID; ++e)
1.275 + s += _lower[e];
1.276 + for (OutEdgeIt e(graph, n); e != INVALID; ++e)
1.277 + s -= _lower[e];
1.278 + supply[n] = s;
1.279 + }
1.280 + valid_supply = true;
1.281 + }
1.282 +
1.283 + /// \brief Simple constructor of the class (without lower bounds).
1.284 + ///
1.285 + /// Simple constructor of the class (without lower bounds).
1.286 + ///
1.287 + /// \param _graph The directed graph the algorithm runs on.
1.288 + /// \param _capacity The capacities (upper bounds) of the edges.
1.289 + /// \param _cost The cost (length) values of the edges.
1.290 + /// \param _s The source node.
1.291 + /// \param _t The target node.
1.292 + /// \param _flow_value The required amount of flow from node \c _s
1.293 + /// to node \c _t (i.e. the supply of \c _s and the demand of
1.294 + /// \c _t).
1.295 + CycleCanceling( const Graph &_graph,
1.296 + const CapacityMap &_capacity,
1.297 + const CostMap &_cost,
1.298 + Node _s, Node _t,
1.299 + Supply _flow_value ) :
1.300 + graph(_graph), lower(NULL), capacity(_capacity), cost(_cost),
1.301 + supply(_graph, 0), demand(supply), flow(_graph, 0),
1.302 + res_graph(_graph, capacity, flow), res_cost(cost)
1.303 + {
1.304 + supply[_s] = _flow_value;
1.305 + supply[_t] = -_flow_value;
1.306 + valid_supply = true;
1.307 + }
1.308 +
1.309 + /// \brief Returns a const reference to the flow map.
1.310 + ///
1.311 + /// Returns a const reference to the flow map.
1.312 + ///
1.313 + /// \pre \ref run() must be called before using this function.
1.314 + const FlowMap& flowMap() const {
1.315 + return flow;
1.316 + }
1.317 +
1.318 + /// \brief Returns the total cost of the found flow.
1.319 + ///
1.320 + /// Returns the total cost of the found flow. The complexity of the
1.321 + /// function is \f$ O(e) \f$.
1.322 + ///
1.323 + /// \pre \ref run() must be called before using this function.
1.324 + Cost totalCost() const {
1.325 + Cost c = 0;
1.326 + for (EdgeIt e(graph); e != INVALID; ++e)
1.327 + c += flow[e] * cost[e];
1.328 + return c;
1.329 + }
1.330 +
1.331 + /// \brief Runs the algorithm.
1.332 + ///
1.333 + /// Runs the algorithm.
1.334 + ///
1.335 + /// \return \c true if a feasible flow can be found.
1.336 + bool run() {
1.337 + return init() && start();
1.338 + }
1.339 +
1.340 + protected:
1.341 +
1.342 + /// \brief Initializes the algorithm.
1.343 + bool init() {
1.344 + // Checking the sum of supply values
1.345 + Supply sum = 0;
1.346 + for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n];
1.347 + if (sum != 0) return false;
1.348 +
1.349 + // Finding a feasible flow
1.350 + Circulation< Graph, Capacity, FlowMap, ConstMap<Edge, Capacity>,
1.351 + CapacityRefMap, DemandMap >
1.352 + circulation( graph, constMap<Edge>((Capacity)0),
1.353 + capacity, demand, flow );
1.354 + return circulation.run() == -1;
1.355 + }
1.356 +
1.357 +#ifdef LIMITED_CYCLE_CANCELING
1.358 + /// \brief Executes a cycle canceling algorithm using
1.359 + /// \ref lemon::BellmanFord "Bellman-Ford" algorithm with limited
1.360 + /// iteration count.
1.361 + bool start() {
1.362 + typename BellmanFord<ResGraph, ResCostMap>::PredMap pred(res_graph);
1.363 + typename ResGraph::template NodeMap<int> visited(res_graph);
1.364 + std::vector<ResEdge> cycle;
1.365 + int node_num = countNodes(graph);
1.366 +
1.367 +#ifdef _DEBUG_ITER_
1.368 + int cycle_num = 0;
1.369 +#endif
1.370 + int length_bound = STARTING_LIMIT;
1.371 + bool optimal = false;
1.372 + while (!optimal) {
1.373 + BellmanFord<ResGraph, ResCostMap> bf(res_graph, res_cost);
1.374 + bf.predMap(pred);
1.375 + bf.init(0);
1.376 + int iter_num = 0;
1.377 + bool cycle_found = false;
1.378 + while (!cycle_found) {
1.379 + int curr_iter_num = iter_num + length_bound <= node_num ?
1.380 + length_bound : node_num - iter_num;
1.381 + iter_num += curr_iter_num;
1.382 + int real_iter_num = curr_iter_num;
1.383 + for (int i = 0; i < curr_iter_num; ++i) {
1.384 + if (bf.processNextWeakRound()) {
1.385 + real_iter_num = i;
1.386 + break;
1.387 + }
1.388 + }
1.389 + if (real_iter_num < curr_iter_num) {
1.390 + optimal = true;
1.391 + break;
1.392 + } else {
1.393 + // Searching for node disjoint negative cycles
1.394 + for (ResNodeIt n(res_graph); n != INVALID; ++n)
1.395 + visited[n] = 0;
1.396 + int id = 0;
1.397 + for (ResNodeIt n(res_graph); n != INVALID; ++n) {
1.398 + if (visited[n] > 0) continue;
1.399 + visited[n] = ++id;
1.400 + ResNode u = pred[n] == INVALID ?
1.401 + INVALID : res_graph.source(pred[n]);
1.402 + while (u != INVALID && visited[u] == 0) {
1.403 + visited[u] = id;
1.404 + u = pred[u] == INVALID ?
1.405 + INVALID : res_graph.source(pred[u]);
1.406 + }
1.407 + if (u != INVALID && visited[u] == id) {
1.408 + // Finding the negative cycle
1.409 + cycle_found = true;
1.410 + cycle.clear();
1.411 + ResEdge e = pred[u];
1.412 + cycle.push_back(e);
1.413 + Capacity d = res_graph.rescap(e);
1.414 + while (res_graph.source(e) != u) {
1.415 + cycle.push_back(e = pred[res_graph.source(e)]);
1.416 + if (res_graph.rescap(e) < d)
1.417 + d = res_graph.rescap(e);
1.418 + }
1.419 +#ifdef _DEBUG_ITER_
1.420 + ++cycle_num;
1.421 +#endif
1.422 + // Augmenting along the cycle
1.423 + for (int i = 0; i < cycle.size(); ++i)
1.424 + res_graph.augment(cycle[i], d);
1.425 +#ifdef _ONLY_ONE_CYCLE_
1.426 + break;
1.427 +#endif
1.428 + }
1.429 + }
1.430 + }
1.431 +
1.432 + if (!cycle_found)
1.433 + length_bound = length_bound * ALPHA_MUL / ALPHA_DIV;
1.434 + }
1.435 + }
1.436 +
1.437 +#ifdef _DEBUG_ITER_
1.438 + std::cout << "Limited cycle canceling algorithm finished. "
1.439 + << "Found " << cycle_num << " negative cycles."
1.440 + << std::endl;
1.441 +#endif
1.442 +
1.443 + // Handling nonzero lower bounds
1.444 + if (lower) {
1.445 + for (EdgeIt e(graph); e != INVALID; ++e)
1.446 + flow[e] += (*lower)[e];
1.447 + }
1.448 + return true;
1.449 + }
1.450 +#endif
1.451 +
1.452 +#ifdef MIN_MEAN_CYCLE_CANCELING
1.453 + /// \brief Executes the minimum mean cycle canceling algorithm
1.454 + /// using \ref lemon::MinMeanCycle "MinMeanCycle" class.
1.455 + bool start() {
1.456 + typedef Path<ResGraph> ResPath;
1.457 + MinMeanCycle<ResGraph, ResCostMap> mmc(res_graph, res_cost);
1.458 + ResPath cycle;
1.459 +
1.460 +#ifdef _DEBUG_ITER_
1.461 + int cycle_num = 0;
1.462 +#endif
1.463 + mmc.cyclePath(cycle).init();
1.464 + if (mmc.findMinMean()) {
1.465 + while (mmc.cycleLength() < 0) {
1.466 +#ifdef _DEBUG_ITER_
1.467 + ++iter;
1.468 +#endif
1.469 + // Finding the cycle
1.470 + mmc.findCycle();
1.471 +
1.472 + // Finding the largest flow amount that can be augmented
1.473 + // along the cycle
1.474 + Capacity delta = 0;
1.475 + for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e) {
1.476 + if (delta == 0 || res_graph.rescap(e) < delta)
1.477 + delta = res_graph.rescap(e);
1.478 + }
1.479 +
1.480 + // Augmenting along the cycle
1.481 + for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e)
1.482 + res_graph.augment(e, delta);
1.483 +
1.484 + // Finding the minimum cycle mean for the modified residual
1.485 + // graph
1.486 + mmc.reset();
1.487 + if (!mmc.findMinMean()) break;
1.488 + }
1.489 + }
1.490 +
1.491 +#ifdef _DEBUG_ITER_
1.492 + std::cout << "Minimum mean cycle canceling algorithm finished. "
1.493 + << "Found " << cycle_num << " negative cycles."
1.494 + << std::endl;
1.495 +#endif
1.496 +
1.497 + // Handling nonzero lower bounds
1.498 + if (lower) {
1.499 + for (EdgeIt e(graph); e != INVALID; ++e)
1.500 + flow[e] += (*lower)[e];
1.501 + }
1.502 + return true;
1.503 + }
1.504 +#endif
1.505 +
1.506 + }; //class CycleCanceling
1.507 +
1.508 + ///@}
1.509 +
1.510 +} //namespace lemon
1.511 +
1.512 +#endif //LEMON_CYCLE_CANCELING_H