deba@2440: /* -*- C++ -*- deba@2440: * deba@2440: * This file is a part of LEMON, a generic C++ optimization library deba@2440: * deba@2440: * Copyright (C) 2003-2007 deba@2440: * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport deba@2440: * (Egervary Research Group on Combinatorial Optimization, EGRES). deba@2440: * deba@2440: * Permission to use, modify and distribute this software is granted deba@2440: * provided that this copyright notice appears in all copies. For deba@2440: * precise terms see the accompanying LICENSE file. deba@2440: * deba@2440: * This software is provided "AS IS" with no warranty of any kind, deba@2440: * express or implied, and with no claim as to its suitability for any deba@2440: * purpose. deba@2440: * deba@2440: */ deba@2440: deba@2440: #ifndef LEMON_CYCLE_CANCELING_H deba@2440: #define LEMON_CYCLE_CANCELING_H deba@2440: deba@2440: /// \ingroup min_cost_flow deba@2440: /// deba@2440: /// \file deba@2440: /// \brief A cycle canceling algorithm for finding a minimum cost flow. deba@2440: deba@2440: #include deba@2440: #include deba@2440: #include deba@2440: deba@2440: /// \brief The used cycle canceling method. deba@2440: #define LIMITED_CYCLE_CANCELING deba@2440: //#define MIN_MEAN_CYCLE_CANCELING deba@2440: deba@2440: #ifdef LIMITED_CYCLE_CANCELING deba@2440: #include deba@2440: /// \brief The maximum number of iterations for the first execution deba@2440: /// of the \ref lemon::BellmanFord "Bellman-Ford" algorithm. deba@2440: #define STARTING_LIMIT 2 deba@2440: /// \brief The iteration limit for the deba@2440: /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by deba@2440: /// ALPHA_MUL % ALPHA_DIV in every round. deba@2440: #define ALPHA_MUL 3 deba@2440: /// \brief The iteration limit for the deba@2440: /// \ref lemon::BellmanFord "Bellman-Ford" algorithm is multiplied by deba@2440: /// ALPHA_MUL % ALPHA_DIV in every round. deba@2440: #define ALPHA_DIV 2 deba@2440: deba@2440: //#define _ONLY_ONE_CYCLE_ deba@2440: //#define _DEBUG_ITER_ deba@2440: #endif deba@2440: deba@2440: #ifdef MIN_MEAN_CYCLE_CANCELING deba@2440: #include deba@2440: #include deba@2440: #endif deba@2440: deba@2440: namespace lemon { deba@2440: deba@2440: /// \addtogroup min_cost_flow deba@2440: /// @{ deba@2440: deba@2440: /// \brief Implementation of a cycle canceling algorithm for finding deba@2440: /// a minimum cost flow. deba@2440: /// deba@2440: /// \ref lemon::CycleCanceling "CycleCanceling" implements a cycle deba@2440: /// canceling algorithm for finding a minimum cost flow. deba@2440: /// deba@2440: /// \param Graph The directed graph type the algorithm runs on. deba@2440: /// \param LowerMap The type of the lower bound map. deba@2440: /// \param CapacityMap The type of the capacity (upper bound) map. deba@2440: /// \param CostMap The type of the cost (length) map. deba@2440: /// \param SupplyMap The type of the supply map. deba@2440: /// deba@2440: /// \warning deba@2440: /// - Edge capacities and costs should be nonnegative integers. deba@2440: /// However \c CostMap::Value should be signed type. deba@2440: /// - Supply values should be integers. deba@2440: /// - \c LowerMap::Value must be convertible to deba@2440: /// \c CapacityMap::Value and \c CapacityMap::Value must be deba@2440: /// convertible to \c SupplyMap::Value. deba@2440: /// deba@2440: /// \author Peter Kovacs deba@2440: deba@2440: template < typename Graph, deba@2440: typename LowerMap = typename Graph::template EdgeMap, deba@2440: typename CapacityMap = LowerMap, deba@2440: typename CostMap = typename Graph::template EdgeMap, deba@2440: typename SupplyMap = typename Graph::template NodeMap deba@2440: > deba@2440: class CycleCanceling deba@2440: { deba@2440: typedef typename Graph::Node Node; deba@2440: typedef typename Graph::NodeIt NodeIt; deba@2440: typedef typename Graph::Edge Edge; deba@2440: typedef typename Graph::EdgeIt EdgeIt; deba@2440: typedef typename Graph::InEdgeIt InEdgeIt; deba@2440: typedef typename Graph::OutEdgeIt OutEdgeIt; deba@2440: deba@2440: typedef typename LowerMap::Value Lower; deba@2440: typedef typename CapacityMap::Value Capacity; deba@2440: typedef typename CostMap::Value Cost; deba@2440: typedef typename SupplyMap::Value Supply; deba@2440: typedef typename Graph::template EdgeMap CapacityRefMap; deba@2440: typedef typename Graph::template NodeMap SupplyRefMap; deba@2440: deba@2440: typedef ResGraphAdaptor< const Graph, Capacity, deba@2440: CapacityRefMap, CapacityRefMap > ResGraph; deba@2440: typedef typename ResGraph::Node ResNode; deba@2440: typedef typename ResGraph::NodeIt ResNodeIt; deba@2440: typedef typename ResGraph::Edge ResEdge; deba@2440: typedef typename ResGraph::EdgeIt ResEdgeIt; deba@2440: deba@2440: public: deba@2440: deba@2440: /// \brief The type of the flow map. deba@2440: typedef CapacityRefMap FlowMap; deba@2440: deba@2440: protected: deba@2440: deba@2440: /// \brief Map adaptor class for demand map. deba@2440: class DemandMap : public MapBase deba@2440: { deba@2440: private: deba@2440: deba@2440: const SupplyMap *map; deba@2440: deba@2440: public: deba@2440: deba@2440: typedef typename MapBase::Value Value; deba@2440: typedef typename MapBase::Key Key; deba@2440: deba@2440: DemandMap(const SupplyMap &_map) : map(&_map) {} deba@2440: deba@2440: Value operator[](const Key &e) const { deba@2440: return -(*map)[e]; deba@2440: } deba@2440: deba@2440: }; //class DemandMap deba@2440: deba@2440: /// \brief Map adaptor class for handling residual edge costs. deba@2440: class ResCostMap : public MapBase deba@2440: { deba@2440: private: deba@2440: deba@2440: const CostMap &cost_map; deba@2440: deba@2440: public: deba@2440: deba@2440: typedef typename MapBase::Value Value; deba@2440: typedef typename MapBase::Key Key; deba@2440: deba@2440: ResCostMap(const CostMap &_cost) : cost_map(_cost) {} deba@2440: deba@2440: Value operator[](const Key &e) const { deba@2440: return ResGraph::forward(e) ? cost_map[e] : -cost_map[e]; deba@2440: } deba@2440: deba@2440: }; //class ResCostMap deba@2440: deba@2440: protected: deba@2440: deba@2440: /// \brief The directed graph the algorithm runs on. deba@2440: const Graph &graph; deba@2440: /// \brief The original lower bound map. deba@2440: const LowerMap *lower; deba@2440: /// \brief The modified capacity map. deba@2440: CapacityRefMap capacity; deba@2440: /// \brief The cost map. deba@2440: const CostMap &cost; deba@2440: /// \brief The modified supply map. deba@2440: SupplyRefMap supply; deba@2440: /// \brief The modified demand map. deba@2440: DemandMap demand; deba@2440: /// \brief The sum of supply values equals zero. deba@2440: bool valid_supply; deba@2440: deba@2440: /// \brief The current flow. deba@2440: FlowMap flow; deba@2440: /// \brief The residual graph. deba@2440: ResGraph res_graph; deba@2440: /// \brief The residual cost map. deba@2440: ResCostMap res_cost; deba@2440: deba@2440: public : deba@2440: deba@2440: /// \brief General constructor of the class (with lower bounds). deba@2440: /// deba@2440: /// General constructor of the class (with lower bounds). deba@2440: /// deba@2440: /// \param _graph The directed graph the algorithm runs on. deba@2440: /// \param _lower The lower bounds of the edges. deba@2440: /// \param _capacity The capacities (upper bounds) of the edges. deba@2440: /// \param _cost The cost (length) values of the edges. deba@2440: /// \param _supply The supply values of the nodes (signed). deba@2440: CycleCanceling( const Graph &_graph, deba@2440: const LowerMap &_lower, deba@2440: const CapacityMap &_capacity, deba@2440: const CostMap &_cost, deba@2440: const SupplyMap &_supply ) : deba@2440: graph(_graph), lower(&_lower), capacity(_graph), cost(_cost), deba@2440: supply(_graph), demand(supply), flow(_graph, 0), deba@2440: res_graph(_graph, capacity, flow), res_cost(cost) deba@2440: { deba@2440: // Removing nonzero lower bounds deba@2440: capacity = subMap(_capacity, _lower); deba@2440: Supply sum = 0; deba@2440: for (NodeIt n(graph); n != INVALID; ++n) { deba@2440: Supply s = _supply[n]; deba@2440: for (InEdgeIt e(graph, n); e != INVALID; ++e) deba@2440: s += _lower[e]; deba@2440: for (OutEdgeIt e(graph, n); e != INVALID; ++e) deba@2440: s -= _lower[e]; deba@2440: sum += (supply[n] = s); deba@2440: } deba@2440: valid_supply = sum == 0; deba@2440: } deba@2440: deba@2440: /// \brief General constructor of the class (without lower bounds). deba@2440: /// deba@2440: /// General constructor of the class (without lower bounds). deba@2440: /// deba@2440: /// \param _graph The directed graph the algorithm runs on. deba@2440: /// \param _capacity The capacities (upper bounds) of the edges. deba@2440: /// \param _cost The cost (length) values of the edges. deba@2440: /// \param _supply The supply values of the nodes (signed). deba@2440: CycleCanceling( const Graph &_graph, deba@2440: const CapacityMap &_capacity, deba@2440: const CostMap &_cost, deba@2440: const SupplyMap &_supply ) : deba@2440: graph(_graph), lower(NULL), capacity(_capacity), cost(_cost), deba@2440: supply(_supply), demand(supply), flow(_graph, 0), deba@2440: res_graph(_graph, capacity, flow), res_cost(cost) deba@2440: { deba@2440: // Checking the sum of supply values deba@2440: Supply sum = 0; deba@2440: for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n]; deba@2440: valid_supply = sum == 0; deba@2440: } deba@2440: deba@2440: deba@2440: /// \brief Simple constructor of the class (with lower bounds). deba@2440: /// deba@2440: /// Simple constructor of the class (with lower bounds). deba@2440: /// deba@2440: /// \param _graph The directed graph the algorithm runs on. deba@2440: /// \param _lower The lower bounds of the edges. deba@2440: /// \param _capacity The capacities (upper bounds) of the edges. deba@2440: /// \param _cost The cost (length) values of the edges. deba@2440: /// \param _s The source node. deba@2440: /// \param _t The target node. deba@2440: /// \param _flow_value The required amount of flow from node \c _s deba@2440: /// to node \c _t (i.e. the supply of \c _s and the demand of deba@2440: /// \c _t). deba@2440: CycleCanceling( const Graph &_graph, deba@2440: const LowerMap &_lower, deba@2440: const CapacityMap &_capacity, deba@2440: const CostMap &_cost, deba@2440: Node _s, Node _t, deba@2440: Supply _flow_value ) : deba@2440: graph(_graph), lower(&_lower), capacity(_graph), cost(_cost), deba@2440: supply(_graph), demand(supply), flow(_graph, 0), deba@2440: res_graph(_graph, capacity, flow), res_cost(cost) deba@2440: { deba@2440: // Removing nonzero lower bounds deba@2440: capacity = subMap(_capacity, _lower); deba@2440: for (NodeIt n(graph); n != INVALID; ++n) { deba@2440: Supply s = 0; deba@2440: if (n == _s) s = _flow_value; deba@2440: if (n == _t) s = -_flow_value; deba@2440: for (InEdgeIt e(graph, n); e != INVALID; ++e) deba@2440: s += _lower[e]; deba@2440: for (OutEdgeIt e(graph, n); e != INVALID; ++e) deba@2440: s -= _lower[e]; deba@2440: supply[n] = s; deba@2440: } deba@2440: valid_supply = true; deba@2440: } deba@2440: deba@2440: /// \brief Simple constructor of the class (without lower bounds). deba@2440: /// deba@2440: /// Simple constructor of the class (without lower bounds). deba@2440: /// deba@2440: /// \param _graph The directed graph the algorithm runs on. deba@2440: /// \param _capacity The capacities (upper bounds) of the edges. deba@2440: /// \param _cost The cost (length) values of the edges. deba@2440: /// \param _s The source node. deba@2440: /// \param _t The target node. deba@2440: /// \param _flow_value The required amount of flow from node \c _s deba@2440: /// to node \c _t (i.e. the supply of \c _s and the demand of deba@2440: /// \c _t). deba@2440: CycleCanceling( const Graph &_graph, deba@2440: const CapacityMap &_capacity, deba@2440: const CostMap &_cost, deba@2440: Node _s, Node _t, deba@2440: Supply _flow_value ) : deba@2440: graph(_graph), lower(NULL), capacity(_capacity), cost(_cost), deba@2440: supply(_graph, 0), demand(supply), flow(_graph, 0), deba@2440: res_graph(_graph, capacity, flow), res_cost(cost) deba@2440: { deba@2440: supply[_s] = _flow_value; deba@2440: supply[_t] = -_flow_value; deba@2440: valid_supply = true; deba@2440: } deba@2440: deba@2440: /// \brief Returns a const reference to the flow map. deba@2440: /// deba@2440: /// Returns a const reference to the flow map. deba@2440: /// deba@2440: /// \pre \ref run() must be called before using this function. deba@2440: const FlowMap& flowMap() const { deba@2440: return flow; deba@2440: } deba@2440: deba@2440: /// \brief Returns the total cost of the found flow. deba@2440: /// deba@2440: /// Returns the total cost of the found flow. The complexity of the deba@2440: /// function is \f$ O(e) \f$. deba@2440: /// deba@2440: /// \pre \ref run() must be called before using this function. deba@2440: Cost totalCost() const { deba@2440: Cost c = 0; deba@2440: for (EdgeIt e(graph); e != INVALID; ++e) deba@2440: c += flow[e] * cost[e]; deba@2440: return c; deba@2440: } deba@2440: deba@2440: /// \brief Runs the algorithm. deba@2440: /// deba@2440: /// Runs the algorithm. deba@2440: /// deba@2440: /// \return \c true if a feasible flow can be found. deba@2440: bool run() { deba@2440: return init() && start(); deba@2440: } deba@2440: deba@2440: protected: deba@2440: deba@2440: /// \brief Initializes the algorithm. deba@2440: bool init() { deba@2440: // Checking the sum of supply values deba@2440: Supply sum = 0; deba@2440: for (NodeIt n(graph); n != INVALID; ++n) sum += supply[n]; deba@2440: if (sum != 0) return false; deba@2440: deba@2440: // Finding a feasible flow deba@2440: Circulation< Graph, Capacity, FlowMap, ConstMap, deba@2440: CapacityRefMap, DemandMap > deba@2440: circulation( graph, constMap((Capacity)0), deba@2440: capacity, demand, flow ); deba@2440: return circulation.run() == -1; deba@2440: } deba@2440: deba@2440: #ifdef LIMITED_CYCLE_CANCELING deba@2440: /// \brief Executes a cycle canceling algorithm using deba@2440: /// \ref lemon::BellmanFord "Bellman-Ford" algorithm with limited deba@2440: /// iteration count. deba@2440: bool start() { deba@2440: typename BellmanFord::PredMap pred(res_graph); deba@2440: typename ResGraph::template NodeMap visited(res_graph); deba@2440: std::vector cycle; deba@2440: int node_num = countNodes(graph); deba@2440: deba@2440: #ifdef _DEBUG_ITER_ deba@2440: int cycle_num = 0; deba@2440: #endif deba@2440: int length_bound = STARTING_LIMIT; deba@2440: bool optimal = false; deba@2440: while (!optimal) { deba@2440: BellmanFord bf(res_graph, res_cost); deba@2440: bf.predMap(pred); deba@2440: bf.init(0); deba@2440: int iter_num = 0; deba@2440: bool cycle_found = false; deba@2440: while (!cycle_found) { deba@2440: int curr_iter_num = iter_num + length_bound <= node_num ? deba@2440: length_bound : node_num - iter_num; deba@2440: iter_num += curr_iter_num; deba@2440: int real_iter_num = curr_iter_num; deba@2440: for (int i = 0; i < curr_iter_num; ++i) { deba@2440: if (bf.processNextWeakRound()) { deba@2440: real_iter_num = i; deba@2440: break; deba@2440: } deba@2440: } deba@2440: if (real_iter_num < curr_iter_num) { deba@2440: optimal = true; deba@2440: break; deba@2440: } else { deba@2440: // Searching for node disjoint negative cycles deba@2440: for (ResNodeIt n(res_graph); n != INVALID; ++n) deba@2440: visited[n] = 0; deba@2440: int id = 0; deba@2440: for (ResNodeIt n(res_graph); n != INVALID; ++n) { deba@2440: if (visited[n] > 0) continue; deba@2440: visited[n] = ++id; deba@2440: ResNode u = pred[n] == INVALID ? deba@2440: INVALID : res_graph.source(pred[n]); deba@2440: while (u != INVALID && visited[u] == 0) { deba@2440: visited[u] = id; deba@2440: u = pred[u] == INVALID ? deba@2440: INVALID : res_graph.source(pred[u]); deba@2440: } deba@2440: if (u != INVALID && visited[u] == id) { deba@2440: // Finding the negative cycle deba@2440: cycle_found = true; deba@2440: cycle.clear(); deba@2440: ResEdge e = pred[u]; deba@2440: cycle.push_back(e); deba@2440: Capacity d = res_graph.rescap(e); deba@2440: while (res_graph.source(e) != u) { deba@2440: cycle.push_back(e = pred[res_graph.source(e)]); deba@2440: if (res_graph.rescap(e) < d) deba@2440: d = res_graph.rescap(e); deba@2440: } deba@2440: #ifdef _DEBUG_ITER_ deba@2440: ++cycle_num; deba@2440: #endif deba@2440: // Augmenting along the cycle deba@2440: for (int i = 0; i < cycle.size(); ++i) deba@2440: res_graph.augment(cycle[i], d); deba@2440: #ifdef _ONLY_ONE_CYCLE_ deba@2440: break; deba@2440: #endif deba@2440: } deba@2440: } deba@2440: } deba@2440: deba@2440: if (!cycle_found) deba@2440: length_bound = length_bound * ALPHA_MUL / ALPHA_DIV; deba@2440: } deba@2440: } deba@2440: deba@2440: #ifdef _DEBUG_ITER_ deba@2440: std::cout << "Limited cycle canceling algorithm finished. " deba@2440: << "Found " << cycle_num << " negative cycles." deba@2440: << std::endl; deba@2440: #endif deba@2440: deba@2440: // Handling nonzero lower bounds deba@2440: if (lower) { deba@2440: for (EdgeIt e(graph); e != INVALID; ++e) deba@2440: flow[e] += (*lower)[e]; deba@2440: } deba@2440: return true; deba@2440: } deba@2440: #endif deba@2440: deba@2440: #ifdef MIN_MEAN_CYCLE_CANCELING deba@2440: /// \brief Executes the minimum mean cycle canceling algorithm deba@2440: /// using \ref lemon::MinMeanCycle "MinMeanCycle" class. deba@2440: bool start() { deba@2440: typedef Path ResPath; deba@2440: MinMeanCycle mmc(res_graph, res_cost); deba@2440: ResPath cycle; deba@2440: deba@2440: #ifdef _DEBUG_ITER_ deba@2440: int cycle_num = 0; deba@2440: #endif deba@2440: mmc.cyclePath(cycle).init(); deba@2440: if (mmc.findMinMean()) { deba@2440: while (mmc.cycleLength() < 0) { deba@2440: #ifdef _DEBUG_ITER_ deba@2440: ++iter; deba@2440: #endif deba@2440: // Finding the cycle deba@2440: mmc.findCycle(); deba@2440: deba@2440: // Finding the largest flow amount that can be augmented deba@2440: // along the cycle deba@2440: Capacity delta = 0; deba@2440: for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e) { deba@2440: if (delta == 0 || res_graph.rescap(e) < delta) deba@2440: delta = res_graph.rescap(e); deba@2440: } deba@2440: deba@2440: // Augmenting along the cycle deba@2440: for (typename ResPath::EdgeIt e(cycle); e != INVALID; ++e) deba@2440: res_graph.augment(e, delta); deba@2440: deba@2440: // Finding the minimum cycle mean for the modified residual deba@2440: // graph deba@2440: mmc.reset(); deba@2440: if (!mmc.findMinMean()) break; deba@2440: } deba@2440: } deba@2440: deba@2440: #ifdef _DEBUG_ITER_ deba@2440: std::cout << "Minimum mean cycle canceling algorithm finished. " deba@2440: << "Found " << cycle_num << " negative cycles." deba@2440: << std::endl; deba@2440: #endif deba@2440: deba@2440: // Handling nonzero lower bounds deba@2440: if (lower) { deba@2440: for (EdgeIt e(graph); e != INVALID; ++e) deba@2440: flow[e] += (*lower)[e]; deba@2440: } deba@2440: return true; deba@2440: } deba@2440: #endif deba@2440: deba@2440: }; //class CycleCanceling deba@2440: deba@2440: ///@} deba@2440: deba@2440: } //namespace lemon deba@2440: deba@2440: #endif //LEMON_CYCLE_CANCELING_H