# HG changeset patch
# User Peter Kovacs <kpeter@inf.elte.hu>
# Date 1258067433 -3600
# Node ID aef153f430e119a38427d6ccc08570cbc0fede3b
# Parent 0643a9c2c3ae4ff91e8faca9061745c21150e0af
Entirely rework cycle canceling algorithms (#180)
- Move the cycle canceling algorithms (CycleCanceling, CancelAndTighten)
into one class (CycleCanceling).
- Add a Method parameter to the run() function to be able to select
the used cycle canceling method.
- Use the new interface similarly to NetworkSimplex.
- Rework the implementations using an efficient internal structure
for handling the residual network.
This improvement made the codes much faster.
- Handle GEQ supply type (LEQ is not supported).
- Handle infinite upper bounds.
- Handle negative costs (for arcs of finite upper bound).
- Extend the documentation.
diff -r 0643a9c2c3ae -r aef153f430e1 lemon/Makefile.am
--- a/lemon/Makefile.am Fri Nov 13 00:09:35 2009 +0100
+++ b/lemon/Makefile.am Fri Nov 13 00:10:33 2009 +0100
@@ -62,7 +62,6 @@
lemon/bin_heap.h \
lemon/binom_heap.h \
lemon/bucket_heap.h \
- lemon/cancel_and_tighten.h \
lemon/capacity_scaling.h \
lemon/cbc.h \
lemon/circulation.h \
diff -r 0643a9c2c3ae -r aef153f430e1 lemon/cancel_and_tighten.h
--- a/lemon/cancel_and_tighten.h Fri Nov 13 00:09:35 2009 +0100
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,797 +0,0 @@
-/* -*- C++ -*-
- *
- * This file is a part of LEMON, a generic C++ optimization library
- *
- * Copyright (C) 2003-2008
- * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
- * (Egervary Research Group on Combinatorial Optimization, EGRES).
- *
- * Permission to use, modify and distribute this software is granted
- * provided that this copyright notice appears in all copies. For
- * precise terms see the accompanying LICENSE file.
- *
- * This software is provided "AS IS" with no warranty of any kind,
- * express or implied, and with no claim as to its suitability for any
- * purpose.
- *
- */
-
-#ifndef LEMON_CANCEL_AND_TIGHTEN_H
-#define LEMON_CANCEL_AND_TIGHTEN_H
-
-/// \ingroup min_cost_flow
-///
-/// \file
-/// \brief Cancel and Tighten algorithm for finding a minimum cost flow.
-
-#include <vector>
-
-#include <lemon/circulation.h>
-#include <lemon/bellman_ford.h>
-#include <lemon/howard.h>
-#include <lemon/adaptors.h>
-#include <lemon/tolerance.h>
-#include <lemon/math.h>
-
-#include <lemon/static_graph.h>
-
-namespace lemon {
-
- /// \addtogroup min_cost_flow
- /// @{
-
- /// \brief Implementation of the Cancel and Tighten algorithm for
- /// finding a minimum cost flow.
- ///
- /// \ref CancelAndTighten implements the Cancel and Tighten algorithm for
- /// finding a minimum cost flow.
- ///
- /// \tparam Digraph The digraph type the algorithm runs on.
- /// \tparam LowerMap The type of the lower bound map.
- /// \tparam CapacityMap The type of the capacity (upper bound) map.
- /// \tparam CostMap The type of the cost (length) map.
- /// \tparam SupplyMap The type of the supply map.
- ///
- /// \warning
- /// - Arc capacities and costs should be \e non-negative \e integers.
- /// - Supply values should be \e signed \e integers.
- /// - The value types of the maps should be convertible to each other.
- /// - \c CostMap::Value must be signed type.
- ///
- /// \author Peter Kovacs
- template < typename Digraph,
- typename LowerMap = typename Digraph::template ArcMap<int>,
- typename CapacityMap = typename Digraph::template ArcMap<int>,
- typename CostMap = typename Digraph::template ArcMap<int>,
- typename SupplyMap = typename Digraph::template NodeMap<int> >
- class CancelAndTighten
- {
- TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
-
- typedef typename CapacityMap::Value Capacity;
- typedef typename CostMap::Value Cost;
- typedef typename SupplyMap::Value Supply;
- typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
- typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
-
- typedef ResidualDigraph< const Digraph,
- CapacityArcMap, CapacityArcMap > ResDigraph;
-
- public:
-
- /// The type of the flow map.
- typedef typename Digraph::template ArcMap<Capacity> FlowMap;
- /// The type of the potential map.
- typedef typename Digraph::template NodeMap<Cost> PotentialMap;
-
- private:
-
- /// \brief Map adaptor class for handling residual arc costs.
- ///
- /// Map adaptor class for handling residual arc costs.
- class ResidualCostMap : public MapBase<typename ResDigraph::Arc, Cost>
- {
- typedef typename ResDigraph::Arc Arc;
-
- private:
-
- const CostMap &_cost_map;
-
- public:
-
- ///\e
- ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
-
- ///\e
- Cost operator[](const Arc &e) const {
- return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
- }
-
- }; //class ResidualCostMap
-
- /// \brief Map adaptor class for handling reduced arc costs.
- ///
- /// Map adaptor class for handling reduced arc costs.
- class ReducedCostMap : public MapBase<Arc, Cost>
- {
- private:
-
- const Digraph &_gr;
- const CostMap &_cost_map;
- const PotentialMap &_pot_map;
-
- public:
-
- ///\e
- ReducedCostMap( const Digraph &gr,
- const CostMap &cost_map,
- const PotentialMap &pot_map ) :
- _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
-
- ///\e
- inline Cost operator[](const Arc &e) const {
- return _cost_map[e] + _pot_map[_gr.source(e)]
- - _pot_map[_gr.target(e)];
- }
-
- }; //class ReducedCostMap
-
- struct BFOperationTraits {
- static double zero() { return 0; }
-
- static double infinity() {
- return std::numeric_limits<double>::infinity();
- }
-
- static double plus(const double& left, const double& right) {
- return left + right;
- }
-
- static bool less(const double& left, const double& right) {
- return left + 1e-6 < right;
- }
- }; // class BFOperationTraits
-
- private:
-
- // The digraph the algorithm runs on
- const Digraph &_graph;
- // The original lower bound map
- const LowerMap *_lower;
- // The modified capacity map
- CapacityArcMap _capacity;
- // The original cost map
- const CostMap &_cost;
- // The modified supply map
- SupplyNodeMap _supply;
- bool _valid_supply;
-
- // Arc map of the current flow
- FlowMap *_flow;
- bool _local_flow;
- // Node map of the current potentials
- PotentialMap *_potential;
- bool _local_potential;
-
- // The residual digraph
- ResDigraph *_res_graph;
- // The residual cost map
- ResidualCostMap _res_cost;
-
- public:
-
- /// \brief General constructor (with lower bounds).
- ///
- /// General constructor (with lower bounds).
- ///
- /// \param digraph The digraph the algorithm runs on.
- /// \param lower The lower bounds of the arcs.
- /// \param capacity The capacities (upper bounds) of the arcs.
- /// \param cost The cost (length) values of the arcs.
- /// \param supply The supply values of the nodes (signed).
- CancelAndTighten( const Digraph &digraph,
- const LowerMap &lower,
- const CapacityMap &capacity,
- const CostMap &cost,
- const SupplyMap &supply ) :
- _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
- _supply(digraph), _flow(NULL), _local_flow(false),
- _potential(NULL), _local_potential(false),
- _res_graph(NULL), _res_cost(_cost)
- {
- // Check the sum of supply values
- Supply sum = 0;
- for (NodeIt n(_graph); n != INVALID; ++n) {
- _supply[n] = supply[n];
- sum += _supply[n];
- }
- _valid_supply = sum == 0;
-
- // Remove non-zero lower bounds
- for (ArcIt e(_graph); e != INVALID; ++e) {
- _capacity[e] = capacity[e];
- if (lower[e] != 0) {
- _capacity[e] -= lower[e];
- _supply[_graph.source(e)] -= lower[e];
- _supply[_graph.target(e)] += lower[e];
- }
- }
- }
-/*
- /// \brief General constructor (without lower bounds).
- ///
- /// General constructor (without lower bounds).
- ///
- /// \param digraph The digraph the algorithm runs on.
- /// \param capacity The capacities (upper bounds) of the arcs.
- /// \param cost The cost (length) values of the arcs.
- /// \param supply The supply values of the nodes (signed).
- CancelAndTighten( const Digraph &digraph,
- const CapacityMap &capacity,
- const CostMap &cost,
- const SupplyMap &supply ) :
- _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
- _supply(supply), _flow(NULL), _local_flow(false),
- _potential(NULL), _local_potential(false),
- _res_graph(NULL), _res_cost(_cost)
- {
- // Check the sum of supply values
- Supply sum = 0;
- for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
- _valid_supply = sum == 0;
- }
-
- /// \brief Simple constructor (with lower bounds).
- ///
- /// Simple constructor (with lower bounds).
- ///
- /// \param digraph The digraph the algorithm runs on.
- /// \param lower The lower bounds of the arcs.
- /// \param capacity The capacities (upper bounds) of the arcs.
- /// \param cost The cost (length) values of the arcs.
- /// \param s The source node.
- /// \param t The target node.
- /// \param flow_value The required amount of flow from node \c s
- /// to node \c t (i.e. the supply of \c s and the demand of \c t).
- CancelAndTighten( const Digraph &digraph,
- const LowerMap &lower,
- const CapacityMap &capacity,
- const CostMap &cost,
- Node s, Node t,
- Supply flow_value ) :
- _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
- _supply(digraph, 0), _flow(NULL), _local_flow(false),
- _potential(NULL), _local_potential(false),
- _res_graph(NULL), _res_cost(_cost)
- {
- // Remove non-zero lower bounds
- _supply[s] = flow_value;
- _supply[t] = -flow_value;
- for (ArcIt e(_graph); e != INVALID; ++e) {
- if (lower[e] != 0) {
- _capacity[e] -= lower[e];
- _supply[_graph.source(e)] -= lower[e];
- _supply[_graph.target(e)] += lower[e];
- }
- }
- _valid_supply = true;
- }
-
- /// \brief Simple constructor (without lower bounds).
- ///
- /// Simple constructor (without lower bounds).
- ///
- /// \param digraph The digraph the algorithm runs on.
- /// \param capacity The capacities (upper bounds) of the arcs.
- /// \param cost The cost (length) values of the arcs.
- /// \param s The source node.
- /// \param t The target node.
- /// \param flow_value The required amount of flow from node \c s
- /// to node \c t (i.e. the supply of \c s and the demand of \c t).
- CancelAndTighten( const Digraph &digraph,
- const CapacityMap &capacity,
- const CostMap &cost,
- Node s, Node t,
- Supply flow_value ) :
- _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
- _supply(digraph, 0), _flow(NULL), _local_flow(false),
- _potential(NULL), _local_potential(false),
- _res_graph(NULL), _res_cost(_cost)
- {
- _supply[s] = flow_value;
- _supply[t] = -flow_value;
- _valid_supply = true;
- }
-*/
- /// Destructor.
- ~CancelAndTighten() {
- if (_local_flow) delete _flow;
- if (_local_potential) delete _potential;
- delete _res_graph;
- }
-
- /// \brief Set the flow map.
- ///
- /// Set the flow map.
- ///
- /// \return \c (*this)
- CancelAndTighten& flowMap(FlowMap &map) {
- if (_local_flow) {
- delete _flow;
- _local_flow = false;
- }
- _flow = ↦
- return *this;
- }
-
- /// \brief Set the potential map.
- ///
- /// Set the potential map.
- ///
- /// \return \c (*this)
- CancelAndTighten& potentialMap(PotentialMap &map) {
- if (_local_potential) {
- delete _potential;
- _local_potential = false;
- }
- _potential = ↦
- return *this;
- }
-
- /// \name Execution control
-
- /// @{
-
- /// \brief Run the algorithm.
- ///
- /// Run the algorithm.
- ///
- /// \return \c true if a feasible flow can be found.
- bool run() {
- return init() && start();
- }
-
- /// @}
-
- /// \name Query Functions
- /// The result of the algorithm can be obtained using these
- /// functions.\n
- /// \ref lemon::CancelAndTighten::run() "run()" must be called before
- /// using them.
-
- /// @{
-
- /// \brief Return a const reference to the arc map storing the
- /// found flow.
- ///
- /// Return a const reference to the arc map storing the found flow.
- ///
- /// \pre \ref run() must be called before using this function.
- const FlowMap& flowMap() const {
- return *_flow;
- }
-
- /// \brief Return a const reference to the node map storing the
- /// found potentials (the dual solution).
- ///
- /// Return a const reference to the node map storing the found
- /// potentials (the dual solution).
- ///
- /// \pre \ref run() must be called before using this function.
- const PotentialMap& potentialMap() const {
- return *_potential;
- }
-
- /// \brief Return the flow on the given arc.
- ///
- /// Return the flow on the given arc.
- ///
- /// \pre \ref run() must be called before using this function.
- Capacity flow(const Arc& arc) const {
- return (*_flow)[arc];
- }
-
- /// \brief Return the potential of the given node.
- ///
- /// Return the potential of the given node.
- ///
- /// \pre \ref run() must be called before using this function.
- Cost potential(const Node& node) const {
- return (*_potential)[node];
- }
-
- /// \brief Return the total cost of the found flow.
- ///
- /// Return the total cost of the found flow. The complexity of the
- /// function is \f$ O(e) \f$.
- ///
- /// \pre \ref run() must be called before using this function.
- Cost totalCost() const {
- Cost c = 0;
- for (ArcIt e(_graph); e != INVALID; ++e)
- c += (*_flow)[e] * _cost[e];
- return c;
- }
-
- /// @}
-
- private:
-
- /// Initialize the algorithm.
- bool init() {
- if (!_valid_supply) return false;
-
- // Initialize flow and potential maps
- if (!_flow) {
- _flow = new FlowMap(_graph);
- _local_flow = true;
- }
- if (!_potential) {
- _potential = new PotentialMap(_graph);
- _local_potential = true;
- }
-
- _res_graph = new ResDigraph(_graph, _capacity, *_flow);
-
- // Find a feasible flow using Circulation
- Circulation< Digraph, ConstMap<Arc, Capacity>,
- CapacityArcMap, SupplyMap >
- circulation( _graph, constMap<Arc>(Capacity(0)),
- _capacity, _supply );
- return circulation.flowMap(*_flow).run();
- }
-
- bool start() {
- const double LIMIT_FACTOR = 0.01;
- const int MIN_LIMIT = 3;
-
- typedef typename Digraph::template NodeMap<double> FloatPotentialMap;
- typedef typename Digraph::template NodeMap<int> LevelMap;
- typedef typename Digraph::template NodeMap<bool> BoolNodeMap;
- typedef typename Digraph::template NodeMap<Node> PredNodeMap;
- typedef typename Digraph::template NodeMap<Arc> PredArcMap;
- typedef typename ResDigraph::template ArcMap<double> ResShiftCostMap;
- FloatPotentialMap pi(_graph);
- LevelMap level(_graph);
- BoolNodeMap reached(_graph);
- BoolNodeMap processed(_graph);
- PredNodeMap pred_node(_graph);
- PredArcMap pred_arc(_graph);
- int node_num = countNodes(_graph);
- typedef std::pair<Arc, bool> pair;
- std::vector<pair> stack(node_num);
- std::vector<Node> proc_vector(node_num);
- ResShiftCostMap shift_cost(*_res_graph);
-
- Tolerance<double> tol;
- tol.epsilon(1e-6);
-
- Timer t1, t2, t3;
- t1.reset();
- t2.reset();
- t3.reset();
-
- // Initialize epsilon and the node potentials
- double epsilon = 0;
- for (ArcIt e(_graph); e != INVALID; ++e) {
- if (_capacity[e] - (*_flow)[e] > 0 && _cost[e] < -epsilon)
- epsilon = -_cost[e];
- else if ((*_flow)[e] > 0 && _cost[e] > epsilon)
- epsilon = _cost[e];
- }
- for (NodeIt v(_graph); v != INVALID; ++v) {
- pi[v] = 0;
- }
-
- // Start phases
- int limit = int(LIMIT_FACTOR * node_num);
- if (limit < MIN_LIMIT) limit = MIN_LIMIT;
- int iter = limit;
- while (epsilon * node_num >= 1) {
- t1.start();
- // Find and cancel cycles in the admissible digraph using DFS
- for (NodeIt n(_graph); n != INVALID; ++n) {
- reached[n] = false;
- processed[n] = false;
- }
- int stack_head = -1;
- int proc_head = -1;
-
- for (NodeIt start(_graph); start != INVALID; ++start) {
- if (reached[start]) continue;
-
- // New start node
- reached[start] = true;
- pred_arc[start] = INVALID;
- pred_node[start] = INVALID;
-
- // Find the first admissible residual outgoing arc
- double p = pi[start];
- Arc e;
- _graph.firstOut(e, start);
- while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
- !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
- _graph.nextOut(e);
- if (e != INVALID) {
- stack[++stack_head] = pair(e, true);
- goto next_step_1;
- }
- _graph.firstIn(e, start);
- while ( e != INVALID && ((*_flow)[e] == 0 ||
- !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
- _graph.nextIn(e);
- if (e != INVALID) {
- stack[++stack_head] = pair(e, false);
- goto next_step_1;
- }
- processed[start] = true;
- proc_vector[++proc_head] = start;
- continue;
- next_step_1:
-
- while (stack_head >= 0) {
- Arc se = stack[stack_head].first;
- bool sf = stack[stack_head].second;
- Node u, v;
- if (sf) {
- u = _graph.source(se);
- v = _graph.target(se);
- } else {
- u = _graph.target(se);
- v = _graph.source(se);
- }
-
- if (!reached[v]) {
- // A new node is reached
- reached[v] = true;
- pred_node[v] = u;
- pred_arc[v] = se;
- // Find the first admissible residual outgoing arc
- double p = pi[v];
- Arc e;
- _graph.firstOut(e, v);
- while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
- !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
- _graph.nextOut(e);
- if (e != INVALID) {
- stack[++stack_head] = pair(e, true);
- goto next_step_2;
- }
- _graph.firstIn(e, v);
- while ( e != INVALID && ((*_flow)[e] == 0 ||
- !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
- _graph.nextIn(e);
- stack[++stack_head] = pair(e, false);
- next_step_2: ;
- } else {
- if (!processed[v]) {
- // A cycle is found
- Node n, w = u;
- Capacity d, delta = sf ? _capacity[se] - (*_flow)[se] :
- (*_flow)[se];
- for (n = u; n != v; n = pred_node[n]) {
- d = _graph.target(pred_arc[n]) == n ?
- _capacity[pred_arc[n]] - (*_flow)[pred_arc[n]] :
- (*_flow)[pred_arc[n]];
- if (d <= delta) {
- delta = d;
- w = pred_node[n];
- }
- }
-
-/*
- std::cout << "CYCLE FOUND: ";
- if (sf)
- std::cout << _cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)];
- else
- std::cout << _graph.id(se) << ":" << -(_cost[se] + pi[_graph.source(se)] - pi[_graph.target(se)]);
- for (n = u; n != v; n = pred_node[n]) {
- if (_graph.target(pred_arc[n]) == n)
- std::cout << " " << _cost[pred_arc[n]] + pi[_graph.source(pred_arc[n])] - pi[_graph.target(pred_arc[n])];
- else
- std::cout << " " << -(_cost[pred_arc[n]] + pi[_graph.source(pred_arc[n])] - pi[_graph.target(pred_arc[n])]);
- }
- std::cout << "\n";
-*/
- // Augment along the cycle
- (*_flow)[se] = sf ? (*_flow)[se] + delta :
- (*_flow)[se] - delta;
- for (n = u; n != v; n = pred_node[n]) {
- if (_graph.target(pred_arc[n]) == n)
- (*_flow)[pred_arc[n]] += delta;
- else
- (*_flow)[pred_arc[n]] -= delta;
- }
- for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
- --stack_head;
- reached[n] = false;
- }
- u = w;
- }
- v = u;
-
- // Find the next admissible residual outgoing arc
- double p = pi[v];
- Arc e = stack[stack_head].first;
- if (!stack[stack_head].second) {
- _graph.nextIn(e);
- goto in_arc_3;
- }
- _graph.nextOut(e);
- while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
- !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
- _graph.nextOut(e);
- if (e != INVALID) {
- stack[stack_head] = pair(e, true);
- goto next_step_3;
- }
- _graph.firstIn(e, v);
- in_arc_3:
- while ( e != INVALID && ((*_flow)[e] == 0 ||
- !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
- _graph.nextIn(e);
- stack[stack_head] = pair(e, false);
- next_step_3: ;
- }
-
- while (stack_head >= 0 && stack[stack_head].first == INVALID) {
- processed[v] = true;
- proc_vector[++proc_head] = v;
- if (--stack_head >= 0) {
- v = stack[stack_head].second ?
- _graph.source(stack[stack_head].first) :
- _graph.target(stack[stack_head].first);
- // Find the next admissible residual outgoing arc
- double p = pi[v];
- Arc e = stack[stack_head].first;
- if (!stack[stack_head].second) {
- _graph.nextIn(e);
- goto in_arc_4;
- }
- _graph.nextOut(e);
- while ( e != INVALID && (_capacity[e] - (*_flow)[e] == 0 ||
- !tol.negative(_cost[e] + p - pi[_graph.target(e)])) )
- _graph.nextOut(e);
- if (e != INVALID) {
- stack[stack_head] = pair(e, true);
- goto next_step_4;
- }
- _graph.firstIn(e, v);
- in_arc_4:
- while ( e != INVALID && ((*_flow)[e] == 0 ||
- !tol.negative(-_cost[e] + p - pi[_graph.source(e)])) )
- _graph.nextIn(e);
- stack[stack_head] = pair(e, false);
- next_step_4: ;
- }
- }
- }
- }
- t1.stop();
-
- // Tighten potentials and epsilon
- if (--iter > 0) {
- // Compute levels
- t2.start();
- for (int i = proc_head; i >= 0; --i) {
- Node v = proc_vector[i];
- double p = pi[v];
- int l = 0;
- for (InArcIt e(_graph, v); e != INVALID; ++e) {
- Node u = _graph.source(e);
- if ( _capacity[e] - (*_flow)[e] > 0 &&
- tol.negative(_cost[e] + pi[u] - p) &&
- level[u] + 1 > l ) l = level[u] + 1;
- }
- for (OutArcIt e(_graph, v); e != INVALID; ++e) {
- Node u = _graph.target(e);
- if ( (*_flow)[e] > 0 &&
- tol.negative(-_cost[e] + pi[u] - p) &&
- level[u] + 1 > l ) l = level[u] + 1;
- }
- level[v] = l;
- }
-
- // Modify potentials
- double p, q = -1;
- for (ArcIt e(_graph); e != INVALID; ++e) {
- Node u = _graph.source(e);
- Node v = _graph.target(e);
- if (_capacity[e] - (*_flow)[e] > 0 && level[u] - level[v] > 0) {
- p = (_cost[e] + pi[u] - pi[v] + epsilon) /
- (level[u] - level[v] + 1);
- if (q < 0 || p < q) q = p;
- }
- else if ((*_flow)[e] > 0 && level[v] - level[u] > 0) {
- p = (-_cost[e] - pi[u] + pi[v] + epsilon) /
- (level[v] - level[u] + 1);
- if (q < 0 || p < q) q = p;
- }
- }
- for (NodeIt v(_graph); v != INVALID; ++v) {
- pi[v] -= q * level[v];
- }
-
- // Modify epsilon
- epsilon = 0;
- for (ArcIt e(_graph); e != INVALID; ++e) {
- double curr = _cost[e] + pi[_graph.source(e)]
- - pi[_graph.target(e)];
- if (_capacity[e] - (*_flow)[e] > 0 && curr < -epsilon)
- epsilon = -curr;
- else if ((*_flow)[e] > 0 && curr > epsilon)
- epsilon = curr;
- }
- t2.stop();
- } else {
- // Set epsilon to the minimum cycle mean
- t3.start();
-
-/**/
- StaticDigraph static_graph;
- typename ResDigraph::template NodeMap<typename StaticDigraph::Node> node_ref(*_res_graph);
- typename ResDigraph::template ArcMap<typename StaticDigraph::Arc> arc_ref(*_res_graph);
- static_graph.build(*_res_graph, node_ref, arc_ref);
- typename StaticDigraph::template NodeMap<double> static_pi(static_graph);
- typename StaticDigraph::template ArcMap<double> static_cost(static_graph);
-
- for (typename ResDigraph::ArcIt e(*_res_graph); e != INVALID; ++e)
- static_cost[arc_ref[e]] = _res_cost[e];
-
- Howard<StaticDigraph, typename StaticDigraph::template ArcMap<double> >
- mmc(static_graph, static_cost);
- mmc.findMinMean();
- epsilon = -mmc.cycleMean();
-/**/
-
-/*
- Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
- mmc.findMinMean();
- epsilon = -mmc.cycleMean();
-*/
-
- // Compute feasible potentials for the current epsilon
- for (typename StaticDigraph::ArcIt e(static_graph); e != INVALID; ++e)
- static_cost[e] += epsilon;
- typename BellmanFord<StaticDigraph, typename StaticDigraph::template ArcMap<double> >::
- template SetDistMap<typename StaticDigraph::template NodeMap<double> >::
- template SetOperationTraits<BFOperationTraits>::Create
- bf(static_graph, static_cost);
- bf.distMap(static_pi).init(0);
- bf.start();
- for (NodeIt n(_graph); n != INVALID; ++n)
- pi[n] = static_pi[node_ref[n]];
-
-/*
- for (typename ResDigraph::ArcIt e(*_res_graph); e != INVALID; ++e)
- shift_cost[e] = _res_cost[e] + epsilon;
- typename BellmanFord<ResDigraph, ResShiftCostMap>::
- template SetDistMap<FloatPotentialMap>::
- template SetOperationTraits<BFOperationTraits>::Create
- bf(*_res_graph, shift_cost);
- bf.distMap(pi).init(0);
- bf.start();
-*/
-
- iter = limit;
- t3.stop();
- }
- }
-
-// std::cout << t1.realTime() << " " << t2.realTime() << " " << t3.realTime() << "\n";
-
- // Handle non-zero lower bounds
- if (_lower) {
- for (ArcIt e(_graph); e != INVALID; ++e)
- (*_flow)[e] += (*_lower)[e];
- }
- return true;
- }
-
- }; //class CancelAndTighten
-
- ///@}
-
-} //namespace lemon
-
-#endif //LEMON_CANCEL_AND_TIGHTEN_H
diff -r 0643a9c2c3ae -r aef153f430e1 lemon/cycle_canceling.h
--- a/lemon/cycle_canceling.h Fri Nov 13 00:09:35 2009 +0100
+++ b/lemon/cycle_canceling.h Fri Nov 13 00:10:33 2009 +0100
@@ -19,441 +19,817 @@
#ifndef LEMON_CYCLE_CANCELING_H
#define LEMON_CYCLE_CANCELING_H
-/// \ingroup min_cost_flow
-///
+/// \ingroup min_cost_flow_algs
/// \file
-/// \brief Cycle-canceling algorithm for finding a minimum cost flow.
+/// \brief Cycle-canceling algorithms for finding a minimum cost flow.
#include <vector>
+#include <limits>
+
+#include <lemon/core.h>
+#include <lemon/maps.h>
+#include <lemon/path.h>
+#include <lemon/math.h>
+#include <lemon/static_graph.h>
#include <lemon/adaptors.h>
-#include <lemon/path.h>
-
#include <lemon/circulation.h>
#include <lemon/bellman_ford.h>
#include <lemon/howard.h>
namespace lemon {
- /// \addtogroup min_cost_flow
+ /// \addtogroup min_cost_flow_algs
/// @{
- /// \brief Implementation of a cycle-canceling algorithm for
- /// finding a minimum cost flow.
+ /// \brief Implementation of cycle-canceling algorithms for
+ /// finding a \ref min_cost_flow "minimum cost flow".
///
- /// \ref CycleCanceling implements a cycle-canceling algorithm for
- /// finding a minimum cost flow.
+ /// \ref CycleCanceling implements three different cycle-canceling
+ /// algorithms for finding a \ref min_cost_flow "minimum cost flow".
+ /// The most efficent one (both theoretically and practically)
+ /// is the \ref CANCEL_AND_TIGHTEN "Cancel and Tighten" algorithm,
+ /// thus it is the default method.
+ /// It is strongly polynomial, but in practice, it is typically much
+ /// slower than the scaling algorithms and NetworkSimplex.
///
- /// \tparam Digraph The digraph type the algorithm runs on.
- /// \tparam LowerMap The type of the lower bound map.
- /// \tparam CapacityMap The type of the capacity (upper bound) map.
- /// \tparam CostMap The type of the cost (length) map.
- /// \tparam SupplyMap The type of the supply map.
+ /// Most of the parameters of the problem (except for the digraph)
+ /// can be given using separate functions, and the algorithm can be
+ /// executed using the \ref run() function. If some parameters are not
+ /// specified, then default values will be used.
///
- /// \warning
- /// - Arc capacities and costs should be \e non-negative \e integers.
- /// - Supply values should be \e signed \e integers.
- /// - The value types of the maps should be convertible to each other.
- /// - \c CostMap::Value must be signed type.
+ /// \tparam GR The digraph type the algorithm runs on.
+ /// \tparam V The number type used for flow amounts, capacity bounds
+ /// and supply values in the algorithm. By default, it is \c int.
+ /// \tparam C The number type used for costs and potentials in the
+ /// algorithm. By default, it is the same as \c V.
///
- /// \note By default the \ref BellmanFord "Bellman-Ford" algorithm is
- /// used for negative cycle detection with limited iteration number.
- /// However \ref CycleCanceling also provides the "Minimum Mean
- /// Cycle-Canceling" algorithm, which is \e strongly \e polynomial,
- /// but rather slower in practice.
- /// To use this version of the algorithm, call \ref run() with \c true
- /// parameter.
+ /// \warning Both number types must be signed and all input data must
+ /// be integer.
+ /// \warning This algorithm does not support negative costs for such
+ /// arcs that have infinite upper bound.
///
- /// \author Peter Kovacs
- template < typename Digraph,
- typename LowerMap = typename Digraph::template ArcMap<int>,
- typename CapacityMap = typename Digraph::template ArcMap<int>,
- typename CostMap = typename Digraph::template ArcMap<int>,
- typename SupplyMap = typename Digraph::template NodeMap<int> >
+ /// \note For more information about the three available methods,
+ /// see \ref Method.
+#ifdef DOXYGEN
+ template <typename GR, typename V, typename C>
+#else
+ template <typename GR, typename V = int, typename C = V>
+#endif
class CycleCanceling
{
- TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
+ public:
- typedef typename CapacityMap::Value Capacity;
- typedef typename CostMap::Value Cost;
- typedef typename SupplyMap::Value Supply;
- typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
- typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
-
- typedef ResidualDigraph< const Digraph,
- CapacityArcMap, CapacityArcMap > ResDigraph;
- typedef typename ResDigraph::Node ResNode;
- typedef typename ResDigraph::NodeIt ResNodeIt;
- typedef typename ResDigraph::Arc ResArc;
- typedef typename ResDigraph::ArcIt ResArcIt;
+ /// The type of the digraph
+ typedef GR Digraph;
+ /// The type of the flow amounts, capacity bounds and supply values
+ typedef V Value;
+ /// The type of the arc costs
+ typedef C Cost;
public:
- /// The type of the flow map.
- typedef typename Digraph::template ArcMap<Capacity> FlowMap;
- /// The type of the potential map.
- typedef typename Digraph::template NodeMap<Cost> PotentialMap;
+ /// \brief Problem type constants for the \c run() function.
+ ///
+ /// Enum type containing the problem type constants that can be
+ /// returned by the \ref run() function of the algorithm.
+ enum ProblemType {
+ /// The problem has no feasible solution (flow).
+ INFEASIBLE,
+ /// The problem has optimal solution (i.e. it is feasible and
+ /// bounded), and the algorithm has found optimal flow and node
+ /// potentials (primal and dual solutions).
+ OPTIMAL,
+ /// The digraph contains an arc of negative cost and infinite
+ /// upper bound. It means that the objective function is unbounded
+ /// on that arc, however, note that it could actually be bounded
+ /// over the feasible flows, but this algroithm cannot handle
+ /// these cases.
+ UNBOUNDED
+ };
+
+ /// \brief Constants for selecting the used method.
+ ///
+ /// Enum type containing constants for selecting the used method
+ /// for the \ref run() function.
+ ///
+ /// \ref CycleCanceling provides three different cycle-canceling
+ /// methods. By default, \ref CANCEL_AND_TIGHTEN "Cancel and Tighten"
+ /// is used, which proved to be the most efficient and the most robust
+ /// on various test inputs.
+ /// However, the other methods can be selected using the \ref run()
+ /// function with the proper parameter.
+ enum Method {
+ /// A simple cycle-canceling method, which uses the
+ /// \ref BellmanFord "Bellman-Ford" algorithm with limited iteration
+ /// number for detecting negative cycles in the residual network.
+ SIMPLE_CYCLE_CANCELING,
+ /// The "Minimum Mean Cycle-Canceling" algorithm, which is a
+ /// well-known strongly polynomial method. It improves along a
+ /// \ref min_mean_cycle "minimum mean cycle" in each iteration.
+ /// Its running time complexity is O(n<sup>2</sup>m<sup>3</sup>log(n)).
+ MINIMUM_MEAN_CYCLE_CANCELING,
+ /// The "Cancel And Tighten" algorithm, which can be viewed as an
+ /// improved version of the previous method.
+ /// It is faster both in theory and in practice, its running time
+ /// complexity is O(n<sup>2</sup>m<sup>2</sup>log(n)).
+ CANCEL_AND_TIGHTEN
+ };
private:
- /// \brief Map adaptor class for handling residual arc costs.
- ///
- /// Map adaptor class for handling residual arc costs.
- class ResidualCostMap : public MapBase<ResArc, Cost>
- {
- private:
+ TEMPLATE_DIGRAPH_TYPEDEFS(GR);
+
+ typedef std::vector<int> IntVector;
+ typedef std::vector<char> CharVector;
+ typedef std::vector<double> DoubleVector;
+ typedef std::vector<Value> ValueVector;
+ typedef std::vector<Cost> CostVector;
- const CostMap &_cost_map;
-
+ private:
+
+ template <typename KT, typename VT>
+ class VectorMap {
public:
-
- ///\e
- ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
-
- ///\e
- Cost operator[](const ResArc &e) const {
- return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
+ typedef KT Key;
+ typedef VT Value;
+
+ VectorMap(std::vector<Value>& v) : _v(v) {}
+
+ const Value& operator[](const Key& key) const {
+ return _v[StaticDigraph::id(key)];
}
- }; //class ResidualCostMap
+ Value& operator[](const Key& key) {
+ return _v[StaticDigraph::id(key)];
+ }
+
+ void set(const Key& key, const Value& val) {
+ _v[StaticDigraph::id(key)] = val;
+ }
+
+ private:
+ std::vector<Value>& _v;
+ };
+
+ typedef VectorMap<StaticDigraph::Node, Cost> CostNodeMap;
+ typedef VectorMap<StaticDigraph::Arc, Cost> CostArcMap;
private:
- // The maximum number of iterations for the first execution of the
- // Bellman-Ford algorithm. It should be at least 2.
- static const int BF_FIRST_LIMIT = 2;
- // The iteration limit for the Bellman-Ford algorithm is multiplied
- // by BF_LIMIT_FACTOR/100 in every round.
- static const int BF_LIMIT_FACTOR = 150;
- private:
+ // Data related to the underlying digraph
+ const GR &_graph;
+ int _node_num;
+ int _arc_num;
+ int _res_node_num;
+ int _res_arc_num;
+ int _root;
- // The digraph the algorithm runs on
- const Digraph &_graph;
- // The original lower bound map
- const LowerMap *_lower;
- // The modified capacity map
- CapacityArcMap _capacity;
- // The original cost map
- const CostMap &_cost;
- // The modified supply map
- SupplyNodeMap _supply;
- bool _valid_supply;
+ // Parameters of the problem
+ bool _have_lower;
+ Value _sum_supply;
- // Arc map of the current flow
- FlowMap *_flow;
- bool _local_flow;
- // Node map of the current potentials
- PotentialMap *_potential;
- bool _local_potential;
+ // Data structures for storing the digraph
+ IntNodeMap _node_id;
+ IntArcMap _arc_idf;
+ IntArcMap _arc_idb;
+ IntVector _first_out;
+ CharVector _forward;
+ IntVector _source;
+ IntVector _target;
+ IntVector _reverse;
- // The residual digraph
- ResDigraph *_res_graph;
- // The residual cost map
- ResidualCostMap _res_cost;
+ // Node and arc data
+ ValueVector _lower;
+ ValueVector _upper;
+ CostVector _cost;
+ ValueVector _supply;
+
+ ValueVector _res_cap;
+ CostVector _pi;
+
+ // Data for a StaticDigraph structure
+ typedef std::pair<int, int> IntPair;
+ StaticDigraph _sgr;
+ std::vector<IntPair> _arc_vec;
+ std::vector<Cost> _cost_vec;
+ IntVector _id_vec;
+ CostArcMap _cost_map;
+ CostNodeMap _pi_map;
+
+ public:
+
+ /// \brief Constant for infinite upper bounds (capacities).
+ ///
+ /// Constant for infinite upper bounds (capacities).
+ /// It is \c std::numeric_limits<Value>::infinity() if available,
+ /// \c std::numeric_limits<Value>::max() otherwise.
+ const Value INF;
public:
- /// \brief General constructor (with lower bounds).
+ /// \brief Constructor.
///
- /// General constructor (with lower bounds).
+ /// The constructor of the class.
///
- /// \param digraph The digraph the algorithm runs on.
- /// \param lower The lower bounds of the arcs.
- /// \param capacity The capacities (upper bounds) of the arcs.
- /// \param cost The cost (length) values of the arcs.
- /// \param supply The supply values of the nodes (signed).
- CycleCanceling( const Digraph &digraph,
- const LowerMap &lower,
- const CapacityMap &capacity,
- const CostMap &cost,
- const SupplyMap &supply ) :
- _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
- _supply(digraph), _flow(NULL), _local_flow(false),
- _potential(NULL), _local_potential(false),
- _res_graph(NULL), _res_cost(_cost)
+ /// \param graph The digraph the algorithm runs on.
+ CycleCanceling(const GR& graph) :
+ _graph(graph), _node_id(graph), _arc_idf(graph), _arc_idb(graph),
+ _cost_map(_cost_vec), _pi_map(_pi),
+ INF(std::numeric_limits<Value>::has_infinity ?
+ std::numeric_limits<Value>::infinity() :
+ std::numeric_limits<Value>::max())
{
- // Check the sum of supply values
- Supply sum = 0;
- for (NodeIt n(_graph); n != INVALID; ++n) {
- _supply[n] = supply[n];
- sum += _supply[n];
+ // Check the number types
+ LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
+ "The flow type of CycleCanceling must be signed");
+ LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
+ "The cost type of CycleCanceling must be signed");
+
+ // Resize vectors
+ _node_num = countNodes(_graph);
+ _arc_num = countArcs(_graph);
+ _res_node_num = _node_num + 1;
+ _res_arc_num = 2 * (_arc_num + _node_num);
+ _root = _node_num;
+
+ _first_out.resize(_res_node_num + 1);
+ _forward.resize(_res_arc_num);
+ _source.resize(_res_arc_num);
+ _target.resize(_res_arc_num);
+ _reverse.resize(_res_arc_num);
+
+ _lower.resize(_res_arc_num);
+ _upper.resize(_res_arc_num);
+ _cost.resize(_res_arc_num);
+ _supply.resize(_res_node_num);
+
+ _res_cap.resize(_res_arc_num);
+ _pi.resize(_res_node_num);
+
+ _arc_vec.reserve(_res_arc_num);
+ _cost_vec.reserve(_res_arc_num);
+ _id_vec.reserve(_res_arc_num);
+
+ // Copy the graph
+ int i = 0, j = 0, k = 2 * _arc_num + _node_num;
+ for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+ _node_id[n] = i;
}
- _valid_supply = sum == 0;
-
- // Remove non-zero lower bounds
- for (ArcIt e(_graph); e != INVALID; ++e) {
- _capacity[e] = capacity[e];
- if (lower[e] != 0) {
- _capacity[e] -= lower[e];
- _supply[_graph.source(e)] -= lower[e];
- _supply[_graph.target(e)] += lower[e];
+ i = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
+ _first_out[i] = j;
+ for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {
+ _arc_idf[a] = j;
+ _forward[j] = true;
+ _source[j] = i;
+ _target[j] = _node_id[_graph.runningNode(a)];
}
+ for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {
+ _arc_idb[a] = j;
+ _forward[j] = false;
+ _source[j] = i;
+ _target[j] = _node_id[_graph.runningNode(a)];
+ }
+ _forward[j] = false;
+ _source[j] = i;
+ _target[j] = _root;
+ _reverse[j] = k;
+ _forward[k] = true;
+ _source[k] = _root;
+ _target[k] = i;
+ _reverse[k] = j;
+ ++j; ++k;
}
- }
-/*
- /// \brief General constructor (without lower bounds).
- ///
- /// General constructor (without lower bounds).
- ///
- /// \param digraph The digraph the algorithm runs on.
- /// \param capacity The capacities (upper bounds) of the arcs.
- /// \param cost The cost (length) values of the arcs.
- /// \param supply The supply values of the nodes (signed).
- CycleCanceling( const Digraph &digraph,
- const CapacityMap &capacity,
- const CostMap &cost,
- const SupplyMap &supply ) :
- _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
- _supply(supply), _flow(NULL), _local_flow(false),
- _potential(NULL), _local_potential(false), _res_graph(NULL),
- _res_cost(_cost)
- {
- // Check the sum of supply values
- Supply sum = 0;
- for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
- _valid_supply = sum == 0;
+ _first_out[i] = j;
+ _first_out[_res_node_num] = k;
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ int fi = _arc_idf[a];
+ int bi = _arc_idb[a];
+ _reverse[fi] = bi;
+ _reverse[bi] = fi;
+ }
+
+ // Reset parameters
+ reset();
}
- /// \brief Simple constructor (with lower bounds).
+ /// \name Parameters
+ /// The parameters of the algorithm can be specified using these
+ /// functions.
+
+ /// @{
+
+ /// \brief Set the lower bounds on the arcs.
///
- /// Simple constructor (with lower bounds).
+ /// This function sets the lower bounds on the arcs.
+ /// If it is not used before calling \ref run(), the lower bounds
+ /// will be set to zero on all arcs.
///
- /// \param digraph The digraph the algorithm runs on.
- /// \param lower The lower bounds of the arcs.
- /// \param capacity The capacities (upper bounds) of the arcs.
- /// \param cost The cost (length) values of the arcs.
- /// \param s The source node.
- /// \param t The target node.
- /// \param flow_value The required amount of flow from node \c s
- /// to node \c t (i.e. the supply of \c s and the demand of \c t).
- CycleCanceling( const Digraph &digraph,
- const LowerMap &lower,
- const CapacityMap &capacity,
- const CostMap &cost,
- Node s, Node t,
- Supply flow_value ) :
- _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
- _supply(digraph, 0), _flow(NULL), _local_flow(false),
- _potential(NULL), _local_potential(false), _res_graph(NULL),
- _res_cost(_cost)
- {
- // Remove non-zero lower bounds
- _supply[s] = flow_value;
- _supply[t] = -flow_value;
- for (ArcIt e(_graph); e != INVALID; ++e) {
- if (lower[e] != 0) {
- _capacity[e] -= lower[e];
- _supply[_graph.source(e)] -= lower[e];
- _supply[_graph.target(e)] += lower[e];
- }
+ /// \param map An arc map storing the lower bounds.
+ /// Its \c Value type must be convertible to the \c Value type
+ /// of the algorithm.
+ ///
+ /// \return <tt>(*this)</tt>
+ template <typename LowerMap>
+ CycleCanceling& lowerMap(const LowerMap& map) {
+ _have_lower = true;
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ _lower[_arc_idf[a]] = map[a];
+ _lower[_arc_idb[a]] = map[a];
}
- _valid_supply = true;
- }
-
- /// \brief Simple constructor (without lower bounds).
- ///
- /// Simple constructor (without lower bounds).
- ///
- /// \param digraph The digraph the algorithm runs on.
- /// \param capacity The capacities (upper bounds) of the arcs.
- /// \param cost The cost (length) values of the arcs.
- /// \param s The source node.
- /// \param t The target node.
- /// \param flow_value The required amount of flow from node \c s
- /// to node \c t (i.e. the supply of \c s and the demand of \c t).
- CycleCanceling( const Digraph &digraph,
- const CapacityMap &capacity,
- const CostMap &cost,
- Node s, Node t,
- Supply flow_value ) :
- _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
- _supply(digraph, 0), _flow(NULL), _local_flow(false),
- _potential(NULL), _local_potential(false), _res_graph(NULL),
- _res_cost(_cost)
- {
- _supply[s] = flow_value;
- _supply[t] = -flow_value;
- _valid_supply = true;
- }
-*/
- /// Destructor.
- ~CycleCanceling() {
- if (_local_flow) delete _flow;
- if (_local_potential) delete _potential;
- delete _res_graph;
- }
-
- /// \brief Set the flow map.
- ///
- /// Set the flow map.
- ///
- /// \return \c (*this)
- CycleCanceling& flowMap(FlowMap &map) {
- if (_local_flow) {
- delete _flow;
- _local_flow = false;
- }
- _flow = ↦
return *this;
}
- /// \brief Set the potential map.
+ /// \brief Set the upper bounds (capacities) on the arcs.
///
- /// Set the potential map.
+ /// This function sets the upper bounds (capacities) on the arcs.
+ /// If it is not used before calling \ref run(), the upper bounds
+ /// will be set to \ref INF on all arcs (i.e. the flow value will be
+ /// unbounded from above).
///
- /// \return \c (*this)
- CycleCanceling& potentialMap(PotentialMap &map) {
- if (_local_potential) {
- delete _potential;
- _local_potential = false;
+ /// \param map An arc map storing the upper bounds.
+ /// Its \c Value type must be convertible to the \c Value type
+ /// of the algorithm.
+ ///
+ /// \return <tt>(*this)</tt>
+ template<typename UpperMap>
+ CycleCanceling& upperMap(const UpperMap& map) {
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ _upper[_arc_idf[a]] = map[a];
}
- _potential = ↦
return *this;
}
+ /// \brief Set the costs of the arcs.
+ ///
+ /// This function sets the costs of the arcs.
+ /// If it is not used before calling \ref run(), the costs
+ /// will be set to \c 1 on all arcs.
+ ///
+ /// \param map An arc map storing the costs.
+ /// Its \c Value type must be convertible to the \c Cost type
+ /// of the algorithm.
+ ///
+ /// \return <tt>(*this)</tt>
+ template<typename CostMap>
+ CycleCanceling& costMap(const CostMap& map) {
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ _cost[_arc_idf[a]] = map[a];
+ _cost[_arc_idb[a]] = -map[a];
+ }
+ return *this;
+ }
+
+ /// \brief Set the supply values of the nodes.
+ ///
+ /// This function sets the supply values of the nodes.
+ /// If neither this function nor \ref stSupply() is used before
+ /// calling \ref run(), the supply of each node will be set to zero.
+ ///
+ /// \param map A node map storing the supply values.
+ /// Its \c Value type must be convertible to the \c Value type
+ /// of the algorithm.
+ ///
+ /// \return <tt>(*this)</tt>
+ template<typename SupplyMap>
+ CycleCanceling& supplyMap(const SupplyMap& map) {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ _supply[_node_id[n]] = map[n];
+ }
+ return *this;
+ }
+
+ /// \brief Set single source and target nodes and a supply value.
+ ///
+ /// This function sets a single source node and a single target node
+ /// and the required flow value.
+ /// If neither this function nor \ref supplyMap() is used before
+ /// calling \ref run(), the supply of each node will be set to zero.
+ ///
+ /// Using this function has the same effect as using \ref supplyMap()
+ /// with such a map in which \c k is assigned to \c s, \c -k is
+ /// assigned to \c t and all other nodes have zero supply value.
+ ///
+ /// \param s The source node.
+ /// \param t The target node.
+ /// \param k The required amount of flow from node \c s to node \c t
+ /// (i.e. the supply of \c s and the demand of \c t).
+ ///
+ /// \return <tt>(*this)</tt>
+ CycleCanceling& stSupply(const Node& s, const Node& t, Value k) {
+ for (int i = 0; i != _res_node_num; ++i) {
+ _supply[i] = 0;
+ }
+ _supply[_node_id[s]] = k;
+ _supply[_node_id[t]] = -k;
+ return *this;
+ }
+
+ /// @}
+
/// \name Execution control
+ /// The algorithm can be executed using \ref run().
/// @{
/// \brief Run the algorithm.
///
- /// Run the algorithm.
+ /// This function runs the algorithm.
+ /// The paramters can be specified using functions \ref lowerMap(),
+ /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply().
+ /// For example,
+ /// \code
+ /// CycleCanceling<ListDigraph> cc(graph);
+ /// cc.lowerMap(lower).upperMap(upper).costMap(cost)
+ /// .supplyMap(sup).run();
+ /// \endcode
///
- /// \param min_mean_cc Set this parameter to \c true to run the
- /// "Minimum Mean Cycle-Canceling" algorithm, which is strongly
- /// polynomial, but rather slower in practice.
+ /// This function can be called more than once. All the parameters
+ /// that have been given are kept for the next call, unless
+ /// \ref reset() is called, thus only the modified parameters
+ /// have to be set again. See \ref reset() for examples.
+ /// However, the underlying digraph must not be modified after this
+ /// class have been constructed, since it copies and extends the graph.
///
- /// \return \c true if a feasible flow can be found.
- bool run(bool min_mean_cc = false) {
- return init() && start(min_mean_cc);
+ /// \param method The cycle-canceling method that will be used.
+ /// For more information, see \ref Method.
+ ///
+ /// \return \c INFEASIBLE if no feasible flow exists,
+ /// \n \c OPTIMAL if the problem has optimal solution
+ /// (i.e. it is feasible and bounded), and the algorithm has found
+ /// optimal flow and node potentials (primal and dual solutions),
+ /// \n \c UNBOUNDED if the digraph contains an arc of negative cost
+ /// and infinite upper bound. It means that the objective function
+ /// is unbounded on that arc, however, note that it could actually be
+ /// bounded over the feasible flows, but this algroithm cannot handle
+ /// these cases.
+ ///
+ /// \see ProblemType, Method
+ ProblemType run(Method method = CANCEL_AND_TIGHTEN) {
+ ProblemType pt = init();
+ if (pt != OPTIMAL) return pt;
+ start(method);
+ return OPTIMAL;
+ }
+
+ /// \brief Reset all the parameters that have been given before.
+ ///
+ /// This function resets all the paramaters that have been given
+ /// before using functions \ref lowerMap(), \ref upperMap(),
+ /// \ref costMap(), \ref supplyMap(), \ref stSupply().
+ ///
+ /// It is useful for multiple run() calls. If this function is not
+ /// used, all the parameters given before are kept for the next
+ /// \ref run() call.
+ /// However, the underlying digraph must not be modified after this
+ /// class have been constructed, since it copies and extends the graph.
+ ///
+ /// For example,
+ /// \code
+ /// CycleCanceling<ListDigraph> cs(graph);
+ ///
+ /// // First run
+ /// cc.lowerMap(lower).upperMap(upper).costMap(cost)
+ /// .supplyMap(sup).run();
+ ///
+ /// // Run again with modified cost map (reset() is not called,
+ /// // so only the cost map have to be set again)
+ /// cost[e] += 100;
+ /// cc.costMap(cost).run();
+ ///
+ /// // Run again from scratch using reset()
+ /// // (the lower bounds will be set to zero on all arcs)
+ /// cc.reset();
+ /// cc.upperMap(capacity).costMap(cost)
+ /// .supplyMap(sup).run();
+ /// \endcode
+ ///
+ /// \return <tt>(*this)</tt>
+ CycleCanceling& reset() {
+ for (int i = 0; i != _res_node_num; ++i) {
+ _supply[i] = 0;
+ }
+ int limit = _first_out[_root];
+ for (int j = 0; j != limit; ++j) {
+ _lower[j] = 0;
+ _upper[j] = INF;
+ _cost[j] = _forward[j] ? 1 : -1;
+ }
+ for (int j = limit; j != _res_arc_num; ++j) {
+ _lower[j] = 0;
+ _upper[j] = INF;
+ _cost[j] = 0;
+ _cost[_reverse[j]] = 0;
+ }
+ _have_lower = false;
+ return *this;
}
/// @}
/// \name Query Functions
- /// The result of the algorithm can be obtained using these
+ /// The results of the algorithm can be obtained using these
/// functions.\n
- /// \ref lemon::CycleCanceling::run() "run()" must be called before
- /// using them.
+ /// The \ref run() function must be called before using them.
/// @{
- /// \brief Return a const reference to the arc map storing the
- /// found flow.
+ /// \brief Return the total cost of the found flow.
///
- /// Return a const reference to the arc map storing the found flow.
+ /// This function returns the total cost of the found flow.
+ /// Its complexity is O(e).
+ ///
+ /// \note The return type of the function can be specified as a
+ /// template parameter. For example,
+ /// \code
+ /// cc.totalCost<double>();
+ /// \endcode
+ /// It is useful if the total cost cannot be stored in the \c Cost
+ /// type of the algorithm, which is the default return type of the
+ /// function.
///
/// \pre \ref run() must be called before using this function.
- const FlowMap& flowMap() const {
- return *_flow;
+ template <typename Number>
+ Number totalCost() const {
+ Number c = 0;
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ int i = _arc_idb[a];
+ c += static_cast<Number>(_res_cap[i]) *
+ (-static_cast<Number>(_cost[i]));
+ }
+ return c;
}
- /// \brief Return a const reference to the node map storing the
- /// found potentials (the dual solution).
- ///
- /// Return a const reference to the node map storing the found
- /// potentials (the dual solution).
- ///
- /// \pre \ref run() must be called before using this function.
- const PotentialMap& potentialMap() const {
- return *_potential;
+#ifndef DOXYGEN
+ Cost totalCost() const {
+ return totalCost<Cost>();
}
+#endif
/// \brief Return the flow on the given arc.
///
- /// Return the flow on the given arc.
+ /// This function returns the flow on the given arc.
///
/// \pre \ref run() must be called before using this function.
- Capacity flow(const Arc& arc) const {
- return (*_flow)[arc];
+ Value flow(const Arc& a) const {
+ return _res_cap[_arc_idb[a]];
}
- /// \brief Return the potential of the given node.
+ /// \brief Return the flow map (the primal solution).
///
- /// Return the potential of the given node.
+ /// This function copies the flow value on each arc into the given
+ /// map. The \c Value type of the algorithm must be convertible to
+ /// the \c Value type of the map.
///
/// \pre \ref run() must be called before using this function.
- Cost potential(const Node& node) const {
- return (*_potential)[node];
+ template <typename FlowMap>
+ void flowMap(FlowMap &map) const {
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ map.set(a, _res_cap[_arc_idb[a]]);
+ }
}
- /// \brief Return the total cost of the found flow.
+ /// \brief Return the potential (dual value) of the given node.
///
- /// Return the total cost of the found flow. The complexity of the
- /// function is \f$ O(e) \f$.
+ /// This function returns the potential (dual value) of the
+ /// given node.
///
/// \pre \ref run() must be called before using this function.
- Cost totalCost() const {
- Cost c = 0;
- for (ArcIt e(_graph); e != INVALID; ++e)
- c += (*_flow)[e] * _cost[e];
- return c;
+ Cost potential(const Node& n) const {
+ return static_cast<Cost>(_pi[_node_id[n]]);
+ }
+
+ /// \brief Return the potential map (the dual solution).
+ ///
+ /// This function copies the potential (dual value) of each node
+ /// into the given map.
+ /// The \c Cost type of the algorithm must be convertible to the
+ /// \c Value type of the map.
+ ///
+ /// \pre \ref run() must be called before using this function.
+ template <typename PotentialMap>
+ void potentialMap(PotentialMap &map) const {
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ map.set(n, static_cast<Cost>(_pi[_node_id[n]]));
+ }
}
/// @}
private:
- /// Initialize the algorithm.
- bool init() {
- if (!_valid_supply) return false;
+ // Initialize the algorithm
+ ProblemType init() {
+ if (_res_node_num <= 1) return INFEASIBLE;
- // Initializing flow and potential maps
- if (!_flow) {
- _flow = new FlowMap(_graph);
- _local_flow = true;
+ // Check the sum of supply values
+ _sum_supply = 0;
+ for (int i = 0; i != _root; ++i) {
+ _sum_supply += _supply[i];
}
- if (!_potential) {
- _potential = new PotentialMap(_graph);
- _local_potential = true;
+ if (_sum_supply > 0) return INFEASIBLE;
+
+
+ // Initialize vectors
+ for (int i = 0; i != _res_node_num; ++i) {
+ _pi[i] = 0;
+ }
+ ValueVector excess(_supply);
+
+ // Remove infinite upper bounds and check negative arcs
+ const Value MAX = std::numeric_limits<Value>::max();
+ int last_out;
+ if (_have_lower) {
+ for (int i = 0; i != _root; ++i) {
+ last_out = _first_out[i+1];
+ for (int j = _first_out[i]; j != last_out; ++j) {
+ if (_forward[j]) {
+ Value c = _cost[j] < 0 ? _upper[j] : _lower[j];
+ if (c >= MAX) return UNBOUNDED;
+ excess[i] -= c;
+ excess[_target[j]] += c;
+ }
+ }
+ }
+ } else {
+ for (int i = 0; i != _root; ++i) {
+ last_out = _first_out[i+1];
+ for (int j = _first_out[i]; j != last_out; ++j) {
+ if (_forward[j] && _cost[j] < 0) {
+ Value c = _upper[j];
+ if (c >= MAX) return UNBOUNDED;
+ excess[i] -= c;
+ excess[_target[j]] += c;
+ }
+ }
+ }
+ }
+ Value ex, max_cap = 0;
+ for (int i = 0; i != _res_node_num; ++i) {
+ ex = excess[i];
+ if (ex < 0) max_cap -= ex;
+ }
+ for (int j = 0; j != _res_arc_num; ++j) {
+ if (_upper[j] >= MAX) _upper[j] = max_cap;
}
- _res_graph = new ResDigraph(_graph, _capacity, *_flow);
+ // Initialize maps for Circulation and remove non-zero lower bounds
+ ConstMap<Arc, Value> low(0);
+ typedef typename Digraph::template ArcMap<Value> ValueArcMap;
+ typedef typename Digraph::template NodeMap<Value> ValueNodeMap;
+ ValueArcMap cap(_graph), flow(_graph);
+ ValueNodeMap sup(_graph);
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ sup[n] = _supply[_node_id[n]];
+ }
+ if (_have_lower) {
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ int j = _arc_idf[a];
+ Value c = _lower[j];
+ cap[a] = _upper[j] - c;
+ sup[_graph.source(a)] -= c;
+ sup[_graph.target(a)] += c;
+ }
+ } else {
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ cap[a] = _upper[_arc_idf[a]];
+ }
+ }
- // Finding a feasible flow using Circulation
- Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
- SupplyMap >
- circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
- _supply );
- return circulation.flowMap(*_flow).run();
+ // Find a feasible flow using Circulation
+ Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
+ circ(_graph, low, cap, sup);
+ if (!circ.flowMap(flow).run()) return INFEASIBLE;
+
+ // Set residual capacities and handle GEQ supply type
+ if (_sum_supply < 0) {
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ Value fa = flow[a];
+ _res_cap[_arc_idf[a]] = cap[a] - fa;
+ _res_cap[_arc_idb[a]] = fa;
+ sup[_graph.source(a)] -= fa;
+ sup[_graph.target(a)] += fa;
+ }
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ excess[_node_id[n]] = sup[n];
+ }
+ for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
+ int u = _target[a];
+ int ra = _reverse[a];
+ _res_cap[a] = -_sum_supply + 1;
+ _res_cap[ra] = -excess[u];
+ _cost[a] = 0;
+ _cost[ra] = 0;
+ }
+ } else {
+ for (ArcIt a(_graph); a != INVALID; ++a) {
+ Value fa = flow[a];
+ _res_cap[_arc_idf[a]] = cap[a] - fa;
+ _res_cap[_arc_idb[a]] = fa;
+ }
+ for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
+ int ra = _reverse[a];
+ _res_cap[a] = 1;
+ _res_cap[ra] = 0;
+ _cost[a] = 0;
+ _cost[ra] = 0;
+ }
+ }
+
+ return OPTIMAL;
+ }
+
+ // Build a StaticDigraph structure containing the current
+ // residual network
+ void buildResidualNetwork() {
+ _arc_vec.clear();
+ _cost_vec.clear();
+ _id_vec.clear();
+ for (int j = 0; j != _res_arc_num; ++j) {
+ if (_res_cap[j] > 0) {
+ _arc_vec.push_back(IntPair(_source[j], _target[j]));
+ _cost_vec.push_back(_cost[j]);
+ _id_vec.push_back(j);
+ }
+ }
+ _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
}
- bool start(bool min_mean_cc) {
- if (min_mean_cc)
- startMinMean();
- else
- start();
+ // Execute the algorithm and transform the results
+ void start(Method method) {
+ // Execute the algorithm
+ switch (method) {
+ case SIMPLE_CYCLE_CANCELING:
+ startSimpleCycleCanceling();
+ break;
+ case MINIMUM_MEAN_CYCLE_CANCELING:
+ startMinMeanCycleCanceling();
+ break;
+ case CANCEL_AND_TIGHTEN:
+ startCancelAndTighten();
+ break;
+ }
- // Handling non-zero lower bounds
- if (_lower) {
- for (ArcIt e(_graph); e != INVALID; ++e)
- (*_flow)[e] += (*_lower)[e];
+ // Compute node potentials
+ if (method != SIMPLE_CYCLE_CANCELING) {
+ buildResidualNetwork();
+ typename BellmanFord<StaticDigraph, CostArcMap>
+ ::template SetDistMap<CostNodeMap>::Create bf(_sgr, _cost_map);
+ bf.distMap(_pi_map);
+ bf.init(0);
+ bf.start();
}
- return true;
+
+ // Handle non-zero lower bounds
+ if (_have_lower) {
+ int limit = _first_out[_root];
+ for (int j = 0; j != limit; ++j) {
+ if (!_forward[j]) _res_cap[j] += _lower[j];
+ }
+ }
}
- /// \brief Execute the algorithm using \ref BellmanFord.
- ///
- /// Execute the algorithm using the \ref BellmanFord
- /// "Bellman-Ford" algorithm for negative cycle detection with
- /// successively larger limit for the number of iterations.
- void start() {
- typename BellmanFord<ResDigraph, ResidualCostMap>::PredMap pred(*_res_graph);
- typename ResDigraph::template NodeMap<int> visited(*_res_graph);
- std::vector<ResArc> cycle;
- int node_num = countNodes(_graph);
+ // Execute the "Simple Cycle Canceling" method
+ void startSimpleCycleCanceling() {
+ // Constants for computing the iteration limits
+ const int BF_FIRST_LIMIT = 2;
+ const double BF_LIMIT_FACTOR = 1.5;
+
+ typedef VectorMap<StaticDigraph::Arc, Value> FilterMap;
+ typedef FilterArcs<StaticDigraph, FilterMap> ResDigraph;
+ typedef VectorMap<StaticDigraph::Node, StaticDigraph::Arc> PredMap;
+ typedef typename BellmanFord<ResDigraph, CostArcMap>
+ ::template SetDistMap<CostNodeMap>
+ ::template SetPredMap<PredMap>::Create BF;
+
+ // Build the residual network
+ _arc_vec.clear();
+ _cost_vec.clear();
+ for (int j = 0; j != _res_arc_num; ++j) {
+ _arc_vec.push_back(IntPair(_source[j], _target[j]));
+ _cost_vec.push_back(_cost[j]);
+ }
+ _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
+
+ FilterMap filter_map(_res_cap);
+ ResDigraph rgr(_sgr, filter_map);
+ std::vector<int> cycle;
+ std::vector<StaticDigraph::Arc> pred(_res_arc_num);
+ PredMap pred_map(pred);
+ BF bf(rgr, _cost_map);
+ bf.distMap(_pi_map).predMap(pred_map);
int length_bound = BF_FIRST_LIMIT;
bool optimal = false;
while (!optimal) {
- BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
- bf.predMap(pred);
bf.init(0);
int iter_num = 0;
bool cycle_found = false;
while (!cycle_found) {
- int curr_iter_num = iter_num + length_bound <= node_num ?
- length_bound : node_num - iter_num;
+ // Perform some iterations of the Bellman-Ford algorithm
+ int curr_iter_num = iter_num + length_bound <= _node_num ?
+ length_bound : _node_num - iter_num;
iter_num += curr_iter_num;
int real_iter_num = curr_iter_num;
for (int i = 0; i < curr_iter_num; ++i) {
@@ -465,89 +841,290 @@
if (real_iter_num < curr_iter_num) {
// Optimal flow is found
optimal = true;
- // Setting node potentials
- for (NodeIt n(_graph); n != INVALID; ++n)
- (*_potential)[n] = bf.dist(n);
break;
} else {
- // Searching for node disjoint negative cycles
- for (ResNodeIt n(*_res_graph); n != INVALID; ++n)
- visited[n] = 0;
+ // Search for node disjoint negative cycles
+ std::vector<int> state(_res_node_num, 0);
int id = 0;
- for (ResNodeIt n(*_res_graph); n != INVALID; ++n) {
- if (visited[n] > 0) continue;
- visited[n] = ++id;
- ResNode u = pred[n] == INVALID ?
- INVALID : _res_graph->source(pred[n]);
- while (u != INVALID && visited[u] == 0) {
- visited[u] = id;
- u = pred[u] == INVALID ?
- INVALID : _res_graph->source(pred[u]);
+ for (int u = 0; u != _res_node_num; ++u) {
+ if (state[u] != 0) continue;
+ ++id;
+ int v = u;
+ for (; v != -1 && state[v] == 0; v = pred[v] == INVALID ?
+ -1 : rgr.id(rgr.source(pred[v]))) {
+ state[v] = id;
}
- if (u != INVALID && visited[u] == id) {
- // Finding the negative cycle
+ if (v != -1 && state[v] == id) {
+ // A negative cycle is found
cycle_found = true;
cycle.clear();
- ResArc e = pred[u];
- cycle.push_back(e);
- Capacity d = _res_graph->residualCapacity(e);
- while (_res_graph->source(e) != u) {
- cycle.push_back(e = pred[_res_graph->source(e)]);
- if (_res_graph->residualCapacity(e) < d)
- d = _res_graph->residualCapacity(e);
+ StaticDigraph::Arc a = pred[v];
+ Value d, delta = _res_cap[rgr.id(a)];
+ cycle.push_back(rgr.id(a));
+ while (rgr.id(rgr.source(a)) != v) {
+ a = pred_map[rgr.source(a)];
+ d = _res_cap[rgr.id(a)];
+ if (d < delta) delta = d;
+ cycle.push_back(rgr.id(a));
}
- // Augmenting along the cycle
- for (int i = 0; i < int(cycle.size()); ++i)
- _res_graph->augment(cycle[i], d);
+ // Augment along the cycle
+ for (int i = 0; i < int(cycle.size()); ++i) {
+ int j = cycle[i];
+ _res_cap[j] -= delta;
+ _res_cap[_reverse[j]] += delta;
+ }
}
}
}
- if (!cycle_found)
- length_bound = length_bound * BF_LIMIT_FACTOR / 100;
+ // Increase iteration limit if no cycle is found
+ if (!cycle_found) {
+ length_bound = static_cast<int>(length_bound * BF_LIMIT_FACTOR);
+ }
}
}
}
- /// \brief Execute the algorithm using \ref Howard.
- ///
- /// Execute the algorithm using \ref Howard for negative
- /// cycle detection.
- void startMinMean() {
- typedef Path<ResDigraph> ResPath;
- Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
- ResPath cycle;
+ // Execute the "Minimum Mean Cycle Canceling" method
+ void startMinMeanCycleCanceling() {
+ typedef SimplePath<StaticDigraph> SPath;
+ typedef typename SPath::ArcIt SPathArcIt;
+ typedef typename Howard<StaticDigraph, CostArcMap>
+ ::template SetPath<SPath>::Create MMC;
+
+ SPath cycle;
+ MMC mmc(_sgr, _cost_map);
+ mmc.cycle(cycle);
+ buildResidualNetwork();
+ while (mmc.findMinMean() && mmc.cycleLength() < 0) {
+ // Find the cycle
+ mmc.findCycle();
- mmc.cycle(cycle);
- if (mmc.findMinMean()) {
- while (mmc.cycleLength() < 0) {
- // Finding the cycle
- mmc.findCycle();
+ // Compute delta value
+ Value delta = INF;
+ for (SPathArcIt a(cycle); a != INVALID; ++a) {
+ Value d = _res_cap[_id_vec[_sgr.id(a)]];
+ if (d < delta) delta = d;
+ }
- // Finding the largest flow amount that can be augmented
- // along the cycle
- Capacity delta = 0;
- for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) {
- if (delta == 0 || _res_graph->residualCapacity(e) < delta)
- delta = _res_graph->residualCapacity(e);
+ // Augment along the cycle
+ for (SPathArcIt a(cycle); a != INVALID; ++a) {
+ int j = _id_vec[_sgr.id(a)];
+ _res_cap[j] -= delta;
+ _res_cap[_reverse[j]] += delta;
+ }
+
+ // Rebuild the residual network
+ buildResidualNetwork();
+ }
+ }
+
+ // Execute the "Cancel And Tighten" method
+ void startCancelAndTighten() {
+ // Constants for the min mean cycle computations
+ const double LIMIT_FACTOR = 1.0;
+ const int MIN_LIMIT = 5;
+
+ // Contruct auxiliary data vectors
+ DoubleVector pi(_res_node_num, 0.0);
+ IntVector level(_res_node_num);
+ CharVector reached(_res_node_num);
+ CharVector processed(_res_node_num);
+ IntVector pred_node(_res_node_num);
+ IntVector pred_arc(_res_node_num);
+ std::vector<int> stack(_res_node_num);
+ std::vector<int> proc_vector(_res_node_num);
+
+ // Initialize epsilon
+ double epsilon = 0;
+ for (int a = 0; a != _res_arc_num; ++a) {
+ if (_res_cap[a] > 0 && -_cost[a] > epsilon)
+ epsilon = -_cost[a];
+ }
+
+ // Start phases
+ Tolerance<double> tol;
+ tol.epsilon(1e-6);
+ int limit = int(LIMIT_FACTOR * std::sqrt(double(_res_node_num)));
+ if (limit < MIN_LIMIT) limit = MIN_LIMIT;
+ int iter = limit;
+ while (epsilon * _res_node_num >= 1) {
+ // Find and cancel cycles in the admissible network using DFS
+ for (int u = 0; u != _res_node_num; ++u) {
+ reached[u] = false;
+ processed[u] = false;
+ }
+ int stack_head = -1;
+ int proc_head = -1;
+ for (int start = 0; start != _res_node_num; ++start) {
+ if (reached[start]) continue;
+
+ // New start node
+ reached[start] = true;
+ pred_arc[start] = -1;
+ pred_node[start] = -1;
+
+ // Find the first admissible outgoing arc
+ double p = pi[start];
+ int a = _first_out[start];
+ int last_out = _first_out[start+1];
+ for (; a != last_out && (_res_cap[a] == 0 ||
+ !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+ if (a == last_out) {
+ processed[start] = true;
+ proc_vector[++proc_head] = start;
+ continue;
+ }
+ stack[++stack_head] = a;
+
+ while (stack_head >= 0) {
+ int sa = stack[stack_head];
+ int u = _source[sa];
+ int v = _target[sa];
+
+ if (!reached[v]) {
+ // A new node is reached
+ reached[v] = true;
+ pred_node[v] = u;
+ pred_arc[v] = sa;
+ p = pi[v];
+ a = _first_out[v];
+ last_out = _first_out[v+1];
+ for (; a != last_out && (_res_cap[a] == 0 ||
+ !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+ stack[++stack_head] = a == last_out ? -1 : a;
+ } else {
+ if (!processed[v]) {
+ // A cycle is found
+ int n, w = u;
+ Value d, delta = _res_cap[sa];
+ for (n = u; n != v; n = pred_node[n]) {
+ d = _res_cap[pred_arc[n]];
+ if (d <= delta) {
+ delta = d;
+ w = pred_node[n];
+ }
+ }
+
+ // Augment along the cycle
+ _res_cap[sa] -= delta;
+ _res_cap[_reverse[sa]] += delta;
+ for (n = u; n != v; n = pred_node[n]) {
+ int pa = pred_arc[n];
+ _res_cap[pa] -= delta;
+ _res_cap[_reverse[pa]] += delta;
+ }
+ for (n = u; stack_head > 0 && n != w; n = pred_node[n]) {
+ --stack_head;
+ reached[n] = false;
+ }
+ u = w;
+ }
+ v = u;
+
+ // Find the next admissible outgoing arc
+ p = pi[v];
+ a = stack[stack_head] + 1;
+ last_out = _first_out[v+1];
+ for (; a != last_out && (_res_cap[a] == 0 ||
+ !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+ stack[stack_head] = a == last_out ? -1 : a;
+ }
+
+ while (stack_head >= 0 && stack[stack_head] == -1) {
+ processed[v] = true;
+ proc_vector[++proc_head] = v;
+ if (--stack_head >= 0) {
+ // Find the next admissible outgoing arc
+ v = _source[stack[stack_head]];
+ p = pi[v];
+ a = stack[stack_head] + 1;
+ last_out = _first_out[v+1];
+ for (; a != last_out && (_res_cap[a] == 0 ||
+ !tol.negative(_cost[a] + p - pi[_target[a]])); ++a) ;
+ stack[stack_head] = a == last_out ? -1 : a;
+ }
+ }
+ }
+ }
+
+ // Tighten potentials and epsilon
+ if (--iter > 0) {
+ for (int u = 0; u != _res_node_num; ++u) {
+ level[u] = 0;
+ }
+ for (int i = proc_head; i > 0; --i) {
+ int u = proc_vector[i];
+ double p = pi[u];
+ int l = level[u] + 1;
+ int last_out = _first_out[u+1];
+ for (int a = _first_out[u]; a != last_out; ++a) {
+ int v = _target[a];
+ if (_res_cap[a] > 0 && tol.negative(_cost[a] + p - pi[v]) &&
+ l > level[v]) level[v] = l;
+ }
}
- // Augmenting along the cycle
- for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e)
- _res_graph->augment(e, delta);
+ // Modify potentials
+ double q = std::numeric_limits<double>::max();
+ for (int u = 0; u != _res_node_num; ++u) {
+ int lu = level[u];
+ double p, pu = pi[u];
+ int last_out = _first_out[u+1];
+ for (int a = _first_out[u]; a != last_out; ++a) {
+ if (_res_cap[a] == 0) continue;
+ int v = _target[a];
+ int ld = lu - level[v];
+ if (ld > 0) {
+ p = (_cost[a] + pu - pi[v] + epsilon) / (ld + 1);
+ if (p < q) q = p;
+ }
+ }
+ }
+ for (int u = 0; u != _res_node_num; ++u) {
+ pi[u] -= q * level[u];
+ }
- // Finding the minimum cycle mean for the modified residual
- // digraph
- if (!mmc.findMinMean()) break;
+ // Modify epsilon
+ epsilon = 0;
+ for (int u = 0; u != _res_node_num; ++u) {
+ double curr, pu = pi[u];
+ int last_out = _first_out[u+1];
+ for (int a = _first_out[u]; a != last_out; ++a) {
+ if (_res_cap[a] == 0) continue;
+ curr = _cost[a] + pu - pi[_target[a]];
+ if (-curr > epsilon) epsilon = -curr;
+ }
+ }
+ } else {
+ typedef Howard<StaticDigraph, CostArcMap> MMC;
+ typedef typename BellmanFord<StaticDigraph, CostArcMap>
+ ::template SetDistMap<CostNodeMap>::Create BF;
+
+ // Set epsilon to the minimum cycle mean
+ buildResidualNetwork();
+ MMC mmc(_sgr, _cost_map);
+ mmc.findMinMean();
+ epsilon = -mmc.cycleMean();
+ Cost cycle_cost = mmc.cycleLength();
+ int cycle_size = mmc.cycleArcNum();
+
+ // Compute feasible potentials for the current epsilon
+ for (int i = 0; i != int(_cost_vec.size()); ++i) {
+ _cost_vec[i] = cycle_size * _cost_vec[i] - cycle_cost;
+ }
+ BF bf(_sgr, _cost_map);
+ bf.distMap(_pi_map);
+ bf.init(0);
+ bf.start();
+ for (int u = 0; u != _res_node_num; ++u) {
+ pi[u] = static_cast<double>(_pi[u]) / cycle_size;
+ }
+
+ iter = limit;
}
}
-
- // Computing node potentials
- BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
- bf.init(0); bf.start();
- for (NodeIt n(_graph); n != INVALID; ++n)
- (*_potential)[n] = bf.dist(n);
}
}; //class CycleCanceling