lemon/cycle_canceling.h
author Peter Kovacs <kpeter@inf.elte.hu>
Fri, 13 Nov 2009 00:09:35 +0100
changeset 814 0643a9c2c3ae
child 815 aef153f430e1
permissions -rw-r--r--
Port cycle canceling algorithms from SVN -r3524 (#180)
kpeter@814
     1
/* -*- C++ -*-
kpeter@814
     2
 *
kpeter@814
     3
 * This file is a part of LEMON, a generic C++ optimization library
kpeter@814
     4
 *
kpeter@814
     5
 * Copyright (C) 2003-2008
kpeter@814
     6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
kpeter@814
     7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
kpeter@814
     8
 *
kpeter@814
     9
 * Permission to use, modify and distribute this software is granted
kpeter@814
    10
 * provided that this copyright notice appears in all copies. For
kpeter@814
    11
 * precise terms see the accompanying LICENSE file.
kpeter@814
    12
 *
kpeter@814
    13
 * This software is provided "AS IS" with no warranty of any kind,
kpeter@814
    14
 * express or implied, and with no claim as to its suitability for any
kpeter@814
    15
 * purpose.
kpeter@814
    16
 *
kpeter@814
    17
 */
kpeter@814
    18
kpeter@814
    19
#ifndef LEMON_CYCLE_CANCELING_H
kpeter@814
    20
#define LEMON_CYCLE_CANCELING_H
kpeter@814
    21
kpeter@814
    22
/// \ingroup min_cost_flow
kpeter@814
    23
///
kpeter@814
    24
/// \file
kpeter@814
    25
/// \brief Cycle-canceling algorithm for finding a minimum cost flow.
kpeter@814
    26
kpeter@814
    27
#include <vector>
kpeter@814
    28
#include <lemon/adaptors.h>
kpeter@814
    29
#include <lemon/path.h>
kpeter@814
    30
kpeter@814
    31
#include <lemon/circulation.h>
kpeter@814
    32
#include <lemon/bellman_ford.h>
kpeter@814
    33
#include <lemon/howard.h>
kpeter@814
    34
kpeter@814
    35
namespace lemon {
kpeter@814
    36
kpeter@814
    37
  /// \addtogroup min_cost_flow
kpeter@814
    38
  /// @{
kpeter@814
    39
kpeter@814
    40
  /// \brief Implementation of a cycle-canceling algorithm for
kpeter@814
    41
  /// finding a minimum cost flow.
kpeter@814
    42
  ///
kpeter@814
    43
  /// \ref CycleCanceling implements a cycle-canceling algorithm for
kpeter@814
    44
  /// finding a minimum cost flow.
kpeter@814
    45
  ///
kpeter@814
    46
  /// \tparam Digraph The digraph type the algorithm runs on.
kpeter@814
    47
  /// \tparam LowerMap The type of the lower bound map.
kpeter@814
    48
  /// \tparam CapacityMap The type of the capacity (upper bound) map.
kpeter@814
    49
  /// \tparam CostMap The type of the cost (length) map.
kpeter@814
    50
  /// \tparam SupplyMap The type of the supply map.
kpeter@814
    51
  ///
kpeter@814
    52
  /// \warning
kpeter@814
    53
  /// - Arc capacities and costs should be \e non-negative \e integers.
kpeter@814
    54
  /// - Supply values should be \e signed \e integers.
kpeter@814
    55
  /// - The value types of the maps should be convertible to each other.
kpeter@814
    56
  /// - \c CostMap::Value must be signed type.
kpeter@814
    57
  ///
kpeter@814
    58
  /// \note By default the \ref BellmanFord "Bellman-Ford" algorithm is
kpeter@814
    59
  /// used for negative cycle detection with limited iteration number.
kpeter@814
    60
  /// However \ref CycleCanceling also provides the "Minimum Mean
kpeter@814
    61
  /// Cycle-Canceling" algorithm, which is \e strongly \e polynomial,
kpeter@814
    62
  /// but rather slower in practice.
kpeter@814
    63
  /// To use this version of the algorithm, call \ref run() with \c true
kpeter@814
    64
  /// parameter.
kpeter@814
    65
  ///
kpeter@814
    66
  /// \author Peter Kovacs
kpeter@814
    67
  template < typename Digraph,
kpeter@814
    68
             typename LowerMap = typename Digraph::template ArcMap<int>,
kpeter@814
    69
             typename CapacityMap = typename Digraph::template ArcMap<int>,
kpeter@814
    70
             typename CostMap = typename Digraph::template ArcMap<int>,
kpeter@814
    71
             typename SupplyMap = typename Digraph::template NodeMap<int> >
kpeter@814
    72
  class CycleCanceling
kpeter@814
    73
  {
kpeter@814
    74
    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);
kpeter@814
    75
kpeter@814
    76
    typedef typename CapacityMap::Value Capacity;
kpeter@814
    77
    typedef typename CostMap::Value Cost;
kpeter@814
    78
    typedef typename SupplyMap::Value Supply;
kpeter@814
    79
    typedef typename Digraph::template ArcMap<Capacity> CapacityArcMap;
kpeter@814
    80
    typedef typename Digraph::template NodeMap<Supply> SupplyNodeMap;
kpeter@814
    81
kpeter@814
    82
    typedef ResidualDigraph< const Digraph,
kpeter@814
    83
      CapacityArcMap, CapacityArcMap > ResDigraph;
kpeter@814
    84
    typedef typename ResDigraph::Node ResNode;
kpeter@814
    85
    typedef typename ResDigraph::NodeIt ResNodeIt;
kpeter@814
    86
    typedef typename ResDigraph::Arc ResArc;
kpeter@814
    87
    typedef typename ResDigraph::ArcIt ResArcIt;
kpeter@814
    88
kpeter@814
    89
  public:
kpeter@814
    90
kpeter@814
    91
    /// The type of the flow map.
kpeter@814
    92
    typedef typename Digraph::template ArcMap<Capacity> FlowMap;
kpeter@814
    93
    /// The type of the potential map.
kpeter@814
    94
    typedef typename Digraph::template NodeMap<Cost> PotentialMap;
kpeter@814
    95
kpeter@814
    96
  private:
kpeter@814
    97
kpeter@814
    98
    /// \brief Map adaptor class for handling residual arc costs.
kpeter@814
    99
    ///
kpeter@814
   100
    /// Map adaptor class for handling residual arc costs.
kpeter@814
   101
    class ResidualCostMap : public MapBase<ResArc, Cost>
kpeter@814
   102
    {
kpeter@814
   103
    private:
kpeter@814
   104
kpeter@814
   105
      const CostMap &_cost_map;
kpeter@814
   106
kpeter@814
   107
    public:
kpeter@814
   108
kpeter@814
   109
      ///\e
kpeter@814
   110
      ResidualCostMap(const CostMap &cost_map) : _cost_map(cost_map) {}
kpeter@814
   111
kpeter@814
   112
      ///\e
kpeter@814
   113
      Cost operator[](const ResArc &e) const {
kpeter@814
   114
        return ResDigraph::forward(e) ? _cost_map[e] : -_cost_map[e];
kpeter@814
   115
      }
kpeter@814
   116
kpeter@814
   117
    }; //class ResidualCostMap
kpeter@814
   118
kpeter@814
   119
  private:
kpeter@814
   120
kpeter@814
   121
    // The maximum number of iterations for the first execution of the
kpeter@814
   122
    // Bellman-Ford algorithm. It should be at least 2.
kpeter@814
   123
    static const int BF_FIRST_LIMIT  = 2;
kpeter@814
   124
    // The iteration limit for the Bellman-Ford algorithm is multiplied
kpeter@814
   125
    // by BF_LIMIT_FACTOR/100 in every round.
kpeter@814
   126
    static const int BF_LIMIT_FACTOR = 150;
kpeter@814
   127
kpeter@814
   128
  private:
kpeter@814
   129
kpeter@814
   130
    // The digraph the algorithm runs on
kpeter@814
   131
    const Digraph &_graph;
kpeter@814
   132
    // The original lower bound map
kpeter@814
   133
    const LowerMap *_lower;
kpeter@814
   134
    // The modified capacity map
kpeter@814
   135
    CapacityArcMap _capacity;
kpeter@814
   136
    // The original cost map
kpeter@814
   137
    const CostMap &_cost;
kpeter@814
   138
    // The modified supply map
kpeter@814
   139
    SupplyNodeMap _supply;
kpeter@814
   140
    bool _valid_supply;
kpeter@814
   141
kpeter@814
   142
    // Arc map of the current flow
kpeter@814
   143
    FlowMap *_flow;
kpeter@814
   144
    bool _local_flow;
kpeter@814
   145
    // Node map of the current potentials
kpeter@814
   146
    PotentialMap *_potential;
kpeter@814
   147
    bool _local_potential;
kpeter@814
   148
kpeter@814
   149
    // The residual digraph
kpeter@814
   150
    ResDigraph *_res_graph;
kpeter@814
   151
    // The residual cost map
kpeter@814
   152
    ResidualCostMap _res_cost;
kpeter@814
   153
kpeter@814
   154
  public:
kpeter@814
   155
kpeter@814
   156
    /// \brief General constructor (with lower bounds).
kpeter@814
   157
    ///
kpeter@814
   158
    /// General constructor (with lower bounds).
kpeter@814
   159
    ///
kpeter@814
   160
    /// \param digraph The digraph the algorithm runs on.
kpeter@814
   161
    /// \param lower The lower bounds of the arcs.
kpeter@814
   162
    /// \param capacity The capacities (upper bounds) of the arcs.
kpeter@814
   163
    /// \param cost The cost (length) values of the arcs.
kpeter@814
   164
    /// \param supply The supply values of the nodes (signed).
kpeter@814
   165
    CycleCanceling( const Digraph &digraph,
kpeter@814
   166
                    const LowerMap &lower,
kpeter@814
   167
                    const CapacityMap &capacity,
kpeter@814
   168
                    const CostMap &cost,
kpeter@814
   169
                    const SupplyMap &supply ) :
kpeter@814
   170
      _graph(digraph), _lower(&lower), _capacity(digraph), _cost(cost),
kpeter@814
   171
      _supply(digraph), _flow(NULL), _local_flow(false),
kpeter@814
   172
      _potential(NULL), _local_potential(false),
kpeter@814
   173
      _res_graph(NULL), _res_cost(_cost)
kpeter@814
   174
    {
kpeter@814
   175
      // Check the sum of supply values
kpeter@814
   176
      Supply sum = 0;
kpeter@814
   177
      for (NodeIt n(_graph); n != INVALID; ++n) {
kpeter@814
   178
        _supply[n] = supply[n];
kpeter@814
   179
        sum += _supply[n];
kpeter@814
   180
      }
kpeter@814
   181
      _valid_supply = sum == 0;
kpeter@814
   182
kpeter@814
   183
      // Remove non-zero lower bounds
kpeter@814
   184
      for (ArcIt e(_graph); e != INVALID; ++e) {
kpeter@814
   185
        _capacity[e] = capacity[e];
kpeter@814
   186
        if (lower[e] != 0) {
kpeter@814
   187
          _capacity[e] -= lower[e];
kpeter@814
   188
          _supply[_graph.source(e)] -= lower[e];
kpeter@814
   189
          _supply[_graph.target(e)] += lower[e];
kpeter@814
   190
        }
kpeter@814
   191
      }
kpeter@814
   192
    }
kpeter@814
   193
/*
kpeter@814
   194
    /// \brief General constructor (without lower bounds).
kpeter@814
   195
    ///
kpeter@814
   196
    /// General constructor (without lower bounds).
kpeter@814
   197
    ///
kpeter@814
   198
    /// \param digraph The digraph the algorithm runs on.
kpeter@814
   199
    /// \param capacity The capacities (upper bounds) of the arcs.
kpeter@814
   200
    /// \param cost The cost (length) values of the arcs.
kpeter@814
   201
    /// \param supply The supply values of the nodes (signed).
kpeter@814
   202
    CycleCanceling( const Digraph &digraph,
kpeter@814
   203
                    const CapacityMap &capacity,
kpeter@814
   204
                    const CostMap &cost,
kpeter@814
   205
                    const SupplyMap &supply ) :
kpeter@814
   206
      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
kpeter@814
   207
      _supply(supply), _flow(NULL), _local_flow(false),
kpeter@814
   208
      _potential(NULL), _local_potential(false), _res_graph(NULL),
kpeter@814
   209
      _res_cost(_cost)
kpeter@814
   210
    {
kpeter@814
   211
      // Check the sum of supply values
kpeter@814
   212
      Supply sum = 0;
kpeter@814
   213
      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
kpeter@814
   214
      _valid_supply = sum == 0;
kpeter@814
   215
    }
kpeter@814
   216
kpeter@814
   217
    /// \brief Simple constructor (with lower bounds).
kpeter@814
   218
    ///
kpeter@814
   219
    /// Simple constructor (with lower bounds).
kpeter@814
   220
    ///
kpeter@814
   221
    /// \param digraph The digraph the algorithm runs on.
kpeter@814
   222
    /// \param lower The lower bounds of the arcs.
kpeter@814
   223
    /// \param capacity The capacities (upper bounds) of the arcs.
kpeter@814
   224
    /// \param cost The cost (length) values of the arcs.
kpeter@814
   225
    /// \param s The source node.
kpeter@814
   226
    /// \param t The target node.
kpeter@814
   227
    /// \param flow_value The required amount of flow from node \c s
kpeter@814
   228
    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
kpeter@814
   229
    CycleCanceling( const Digraph &digraph,
kpeter@814
   230
                    const LowerMap &lower,
kpeter@814
   231
                    const CapacityMap &capacity,
kpeter@814
   232
                    const CostMap &cost,
kpeter@814
   233
                    Node s, Node t,
kpeter@814
   234
                    Supply flow_value ) :
kpeter@814
   235
      _graph(digraph), _lower(&lower), _capacity(capacity), _cost(cost),
kpeter@814
   236
      _supply(digraph, 0), _flow(NULL), _local_flow(false),
kpeter@814
   237
      _potential(NULL), _local_potential(false), _res_graph(NULL),
kpeter@814
   238
      _res_cost(_cost)
kpeter@814
   239
    {
kpeter@814
   240
      // Remove non-zero lower bounds
kpeter@814
   241
      _supply[s] =  flow_value;
kpeter@814
   242
      _supply[t] = -flow_value;
kpeter@814
   243
      for (ArcIt e(_graph); e != INVALID; ++e) {
kpeter@814
   244
        if (lower[e] != 0) {
kpeter@814
   245
          _capacity[e] -= lower[e];
kpeter@814
   246
          _supply[_graph.source(e)] -= lower[e];
kpeter@814
   247
          _supply[_graph.target(e)] += lower[e];
kpeter@814
   248
        }
kpeter@814
   249
      }
kpeter@814
   250
      _valid_supply = true;
kpeter@814
   251
    }
kpeter@814
   252
kpeter@814
   253
    /// \brief Simple constructor (without lower bounds).
kpeter@814
   254
    ///
kpeter@814
   255
    /// Simple constructor (without lower bounds).
kpeter@814
   256
    ///
kpeter@814
   257
    /// \param digraph The digraph the algorithm runs on.
kpeter@814
   258
    /// \param capacity The capacities (upper bounds) of the arcs.
kpeter@814
   259
    /// \param cost The cost (length) values of the arcs.
kpeter@814
   260
    /// \param s The source node.
kpeter@814
   261
    /// \param t The target node.
kpeter@814
   262
    /// \param flow_value The required amount of flow from node \c s
kpeter@814
   263
    /// to node \c t (i.e. the supply of \c s and the demand of \c t).
kpeter@814
   264
    CycleCanceling( const Digraph &digraph,
kpeter@814
   265
                    const CapacityMap &capacity,
kpeter@814
   266
                    const CostMap &cost,
kpeter@814
   267
                    Node s, Node t,
kpeter@814
   268
                    Supply flow_value ) :
kpeter@814
   269
      _graph(digraph), _lower(NULL), _capacity(capacity), _cost(cost),
kpeter@814
   270
      _supply(digraph, 0), _flow(NULL), _local_flow(false),
kpeter@814
   271
      _potential(NULL), _local_potential(false), _res_graph(NULL),
kpeter@814
   272
      _res_cost(_cost)
kpeter@814
   273
    {
kpeter@814
   274
      _supply[s] =  flow_value;
kpeter@814
   275
      _supply[t] = -flow_value;
kpeter@814
   276
      _valid_supply = true;
kpeter@814
   277
    }
kpeter@814
   278
*/
kpeter@814
   279
    /// Destructor.
kpeter@814
   280
    ~CycleCanceling() {
kpeter@814
   281
      if (_local_flow) delete _flow;
kpeter@814
   282
      if (_local_potential) delete _potential;
kpeter@814
   283
      delete _res_graph;
kpeter@814
   284
    }
kpeter@814
   285
kpeter@814
   286
    /// \brief Set the flow map.
kpeter@814
   287
    ///
kpeter@814
   288
    /// Set the flow map.
kpeter@814
   289
    ///
kpeter@814
   290
    /// \return \c (*this)
kpeter@814
   291
    CycleCanceling& flowMap(FlowMap &map) {
kpeter@814
   292
      if (_local_flow) {
kpeter@814
   293
        delete _flow;
kpeter@814
   294
        _local_flow = false;
kpeter@814
   295
      }
kpeter@814
   296
      _flow = &map;
kpeter@814
   297
      return *this;
kpeter@814
   298
    }
kpeter@814
   299
kpeter@814
   300
    /// \brief Set the potential map.
kpeter@814
   301
    ///
kpeter@814
   302
    /// Set the potential map.
kpeter@814
   303
    ///
kpeter@814
   304
    /// \return \c (*this)
kpeter@814
   305
    CycleCanceling& potentialMap(PotentialMap &map) {
kpeter@814
   306
      if (_local_potential) {
kpeter@814
   307
        delete _potential;
kpeter@814
   308
        _local_potential = false;
kpeter@814
   309
      }
kpeter@814
   310
      _potential = &map;
kpeter@814
   311
      return *this;
kpeter@814
   312
    }
kpeter@814
   313
kpeter@814
   314
    /// \name Execution control
kpeter@814
   315
kpeter@814
   316
    /// @{
kpeter@814
   317
kpeter@814
   318
    /// \brief Run the algorithm.
kpeter@814
   319
    ///
kpeter@814
   320
    /// Run the algorithm.
kpeter@814
   321
    ///
kpeter@814
   322
    /// \param min_mean_cc Set this parameter to \c true to run the
kpeter@814
   323
    /// "Minimum Mean Cycle-Canceling" algorithm, which is strongly
kpeter@814
   324
    /// polynomial, but rather slower in practice.
kpeter@814
   325
    ///
kpeter@814
   326
    /// \return \c true if a feasible flow can be found.
kpeter@814
   327
    bool run(bool min_mean_cc = false) {
kpeter@814
   328
      return init() && start(min_mean_cc);
kpeter@814
   329
    }
kpeter@814
   330
kpeter@814
   331
    /// @}
kpeter@814
   332
kpeter@814
   333
    /// \name Query Functions
kpeter@814
   334
    /// The result of the algorithm can be obtained using these
kpeter@814
   335
    /// functions.\n
kpeter@814
   336
    /// \ref lemon::CycleCanceling::run() "run()" must be called before
kpeter@814
   337
    /// using them.
kpeter@814
   338
kpeter@814
   339
    /// @{
kpeter@814
   340
kpeter@814
   341
    /// \brief Return a const reference to the arc map storing the
kpeter@814
   342
    /// found flow.
kpeter@814
   343
    ///
kpeter@814
   344
    /// Return a const reference to the arc map storing the found flow.
kpeter@814
   345
    ///
kpeter@814
   346
    /// \pre \ref run() must be called before using this function.
kpeter@814
   347
    const FlowMap& flowMap() const {
kpeter@814
   348
      return *_flow;
kpeter@814
   349
    }
kpeter@814
   350
kpeter@814
   351
    /// \brief Return a const reference to the node map storing the
kpeter@814
   352
    /// found potentials (the dual solution).
kpeter@814
   353
    ///
kpeter@814
   354
    /// Return a const reference to the node map storing the found
kpeter@814
   355
    /// potentials (the dual solution).
kpeter@814
   356
    ///
kpeter@814
   357
    /// \pre \ref run() must be called before using this function.
kpeter@814
   358
    const PotentialMap& potentialMap() const {
kpeter@814
   359
      return *_potential;
kpeter@814
   360
    }
kpeter@814
   361
kpeter@814
   362
    /// \brief Return the flow on the given arc.
kpeter@814
   363
    ///
kpeter@814
   364
    /// Return the flow on the given arc.
kpeter@814
   365
    ///
kpeter@814
   366
    /// \pre \ref run() must be called before using this function.
kpeter@814
   367
    Capacity flow(const Arc& arc) const {
kpeter@814
   368
      return (*_flow)[arc];
kpeter@814
   369
    }
kpeter@814
   370
kpeter@814
   371
    /// \brief Return the potential of the given node.
kpeter@814
   372
    ///
kpeter@814
   373
    /// Return the potential of the given node.
kpeter@814
   374
    ///
kpeter@814
   375
    /// \pre \ref run() must be called before using this function.
kpeter@814
   376
    Cost potential(const Node& node) const {
kpeter@814
   377
      return (*_potential)[node];
kpeter@814
   378
    }
kpeter@814
   379
kpeter@814
   380
    /// \brief Return the total cost of the found flow.
kpeter@814
   381
    ///
kpeter@814
   382
    /// Return the total cost of the found flow. The complexity of the
kpeter@814
   383
    /// function is \f$ O(e) \f$.
kpeter@814
   384
    ///
kpeter@814
   385
    /// \pre \ref run() must be called before using this function.
kpeter@814
   386
    Cost totalCost() const {
kpeter@814
   387
      Cost c = 0;
kpeter@814
   388
      for (ArcIt e(_graph); e != INVALID; ++e)
kpeter@814
   389
        c += (*_flow)[e] * _cost[e];
kpeter@814
   390
      return c;
kpeter@814
   391
    }
kpeter@814
   392
kpeter@814
   393
    /// @}
kpeter@814
   394
kpeter@814
   395
  private:
kpeter@814
   396
kpeter@814
   397
    /// Initialize the algorithm.
kpeter@814
   398
    bool init() {
kpeter@814
   399
      if (!_valid_supply) return false;
kpeter@814
   400
kpeter@814
   401
      // Initializing flow and potential maps
kpeter@814
   402
      if (!_flow) {
kpeter@814
   403
        _flow = new FlowMap(_graph);
kpeter@814
   404
        _local_flow = true;
kpeter@814
   405
      }
kpeter@814
   406
      if (!_potential) {
kpeter@814
   407
        _potential = new PotentialMap(_graph);
kpeter@814
   408
        _local_potential = true;
kpeter@814
   409
      }
kpeter@814
   410
kpeter@814
   411
      _res_graph = new ResDigraph(_graph, _capacity, *_flow);
kpeter@814
   412
kpeter@814
   413
      // Finding a feasible flow using Circulation
kpeter@814
   414
      Circulation< Digraph, ConstMap<Arc, Capacity>, CapacityArcMap,
kpeter@814
   415
                   SupplyMap >
kpeter@814
   416
        circulation( _graph, constMap<Arc>(Capacity(0)), _capacity,
kpeter@814
   417
                     _supply );
kpeter@814
   418
      return circulation.flowMap(*_flow).run();
kpeter@814
   419
    }
kpeter@814
   420
kpeter@814
   421
    bool start(bool min_mean_cc) {
kpeter@814
   422
      if (min_mean_cc)
kpeter@814
   423
        startMinMean();
kpeter@814
   424
      else
kpeter@814
   425
        start();
kpeter@814
   426
kpeter@814
   427
      // Handling non-zero lower bounds
kpeter@814
   428
      if (_lower) {
kpeter@814
   429
        for (ArcIt e(_graph); e != INVALID; ++e)
kpeter@814
   430
          (*_flow)[e] += (*_lower)[e];
kpeter@814
   431
      }
kpeter@814
   432
      return true;
kpeter@814
   433
    }
kpeter@814
   434
kpeter@814
   435
    /// \brief Execute the algorithm using \ref BellmanFord.
kpeter@814
   436
    ///
kpeter@814
   437
    /// Execute the algorithm using the \ref BellmanFord
kpeter@814
   438
    /// "Bellman-Ford" algorithm for negative cycle detection with
kpeter@814
   439
    /// successively larger limit for the number of iterations.
kpeter@814
   440
    void start() {
kpeter@814
   441
      typename BellmanFord<ResDigraph, ResidualCostMap>::PredMap pred(*_res_graph);
kpeter@814
   442
      typename ResDigraph::template NodeMap<int> visited(*_res_graph);
kpeter@814
   443
      std::vector<ResArc> cycle;
kpeter@814
   444
      int node_num = countNodes(_graph);
kpeter@814
   445
kpeter@814
   446
      int length_bound = BF_FIRST_LIMIT;
kpeter@814
   447
      bool optimal = false;
kpeter@814
   448
      while (!optimal) {
kpeter@814
   449
        BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
kpeter@814
   450
        bf.predMap(pred);
kpeter@814
   451
        bf.init(0);
kpeter@814
   452
        int iter_num = 0;
kpeter@814
   453
        bool cycle_found = false;
kpeter@814
   454
        while (!cycle_found) {
kpeter@814
   455
          int curr_iter_num = iter_num + length_bound <= node_num ?
kpeter@814
   456
                              length_bound : node_num - iter_num;
kpeter@814
   457
          iter_num += curr_iter_num;
kpeter@814
   458
          int real_iter_num = curr_iter_num;
kpeter@814
   459
          for (int i = 0; i < curr_iter_num; ++i) {
kpeter@814
   460
            if (bf.processNextWeakRound()) {
kpeter@814
   461
              real_iter_num = i;
kpeter@814
   462
              break;
kpeter@814
   463
            }
kpeter@814
   464
          }
kpeter@814
   465
          if (real_iter_num < curr_iter_num) {
kpeter@814
   466
            // Optimal flow is found
kpeter@814
   467
            optimal = true;
kpeter@814
   468
            // Setting node potentials
kpeter@814
   469
            for (NodeIt n(_graph); n != INVALID; ++n)
kpeter@814
   470
              (*_potential)[n] = bf.dist(n);
kpeter@814
   471
            break;
kpeter@814
   472
          } else {
kpeter@814
   473
            // Searching for node disjoint negative cycles
kpeter@814
   474
            for (ResNodeIt n(*_res_graph); n != INVALID; ++n)
kpeter@814
   475
              visited[n] = 0;
kpeter@814
   476
            int id = 0;
kpeter@814
   477
            for (ResNodeIt n(*_res_graph); n != INVALID; ++n) {
kpeter@814
   478
              if (visited[n] > 0) continue;
kpeter@814
   479
              visited[n] = ++id;
kpeter@814
   480
              ResNode u = pred[n] == INVALID ?
kpeter@814
   481
                          INVALID : _res_graph->source(pred[n]);
kpeter@814
   482
              while (u != INVALID && visited[u] == 0) {
kpeter@814
   483
                visited[u] = id;
kpeter@814
   484
                u = pred[u] == INVALID ?
kpeter@814
   485
                    INVALID : _res_graph->source(pred[u]);
kpeter@814
   486
              }
kpeter@814
   487
              if (u != INVALID && visited[u] == id) {
kpeter@814
   488
                // Finding the negative cycle
kpeter@814
   489
                cycle_found = true;
kpeter@814
   490
                cycle.clear();
kpeter@814
   491
                ResArc e = pred[u];
kpeter@814
   492
                cycle.push_back(e);
kpeter@814
   493
                Capacity d = _res_graph->residualCapacity(e);
kpeter@814
   494
                while (_res_graph->source(e) != u) {
kpeter@814
   495
                  cycle.push_back(e = pred[_res_graph->source(e)]);
kpeter@814
   496
                  if (_res_graph->residualCapacity(e) < d)
kpeter@814
   497
                    d = _res_graph->residualCapacity(e);
kpeter@814
   498
                }
kpeter@814
   499
kpeter@814
   500
                // Augmenting along the cycle
kpeter@814
   501
                for (int i = 0; i < int(cycle.size()); ++i)
kpeter@814
   502
                  _res_graph->augment(cycle[i], d);
kpeter@814
   503
              }
kpeter@814
   504
            }
kpeter@814
   505
          }
kpeter@814
   506
kpeter@814
   507
          if (!cycle_found)
kpeter@814
   508
            length_bound = length_bound * BF_LIMIT_FACTOR / 100;
kpeter@814
   509
        }
kpeter@814
   510
      }
kpeter@814
   511
    }
kpeter@814
   512
kpeter@814
   513
    /// \brief Execute the algorithm using \ref Howard.
kpeter@814
   514
    ///
kpeter@814
   515
    /// Execute the algorithm using \ref Howard for negative
kpeter@814
   516
    /// cycle detection.
kpeter@814
   517
    void startMinMean() {
kpeter@814
   518
      typedef Path<ResDigraph> ResPath;
kpeter@814
   519
      Howard<ResDigraph, ResidualCostMap> mmc(*_res_graph, _res_cost);
kpeter@814
   520
      ResPath cycle;
kpeter@814
   521
kpeter@814
   522
      mmc.cycle(cycle);
kpeter@814
   523
      if (mmc.findMinMean()) {
kpeter@814
   524
        while (mmc.cycleLength() < 0) {
kpeter@814
   525
          // Finding the cycle
kpeter@814
   526
          mmc.findCycle();
kpeter@814
   527
kpeter@814
   528
          // Finding the largest flow amount that can be augmented
kpeter@814
   529
          // along the cycle
kpeter@814
   530
          Capacity delta = 0;
kpeter@814
   531
          for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e) {
kpeter@814
   532
            if (delta == 0 || _res_graph->residualCapacity(e) < delta)
kpeter@814
   533
              delta = _res_graph->residualCapacity(e);
kpeter@814
   534
          }
kpeter@814
   535
kpeter@814
   536
          // Augmenting along the cycle
kpeter@814
   537
          for (typename ResPath::ArcIt e(cycle); e != INVALID; ++e)
kpeter@814
   538
            _res_graph->augment(e, delta);
kpeter@814
   539
kpeter@814
   540
          // Finding the minimum cycle mean for the modified residual
kpeter@814
   541
          // digraph
kpeter@814
   542
          if (!mmc.findMinMean()) break;
kpeter@814
   543
        }
kpeter@814
   544
      }
kpeter@814
   545
kpeter@814
   546
      // Computing node potentials
kpeter@814
   547
      BellmanFord<ResDigraph, ResidualCostMap> bf(*_res_graph, _res_cost);
kpeter@814
   548
      bf.init(0); bf.start();
kpeter@814
   549
      for (NodeIt n(_graph); n != INVALID; ++n)
kpeter@814
   550
        (*_potential)[n] = bf.dist(n);
kpeter@814
   551
    }
kpeter@814
   552
kpeter@814
   553
  }; //class CycleCanceling
kpeter@814
   554
kpeter@814
   555
  ///@}
kpeter@814
   556
kpeter@814
   557
} //namespace lemon
kpeter@814
   558
kpeter@814
   559
#endif //LEMON_CYCLE_CANCELING_H