COIN-OR::LEMON - Graph Library

Changeset 2634:e98bbe64cca4 in lemon-0.x


Ignore:
Timestamp:
02/06/09 22:52:34 (11 years ago)
Author:
Peter Kovacs
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3519
Message:

Rework Network Simplex
Use simpler and faster graph implementation instead of SmartGraph?

File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/network_simplex.h

    r2630 r2634  
    2929#include <algorithm>
    3030
    31 #include <lemon/graph_adaptor.h>
    3231#include <lemon/graph_utils.h>
    33 #include <lemon/smart_graph.h>
    3432#include <lemon/math.h>
    3533
     
    7371  class NetworkSimplex
    7472  {
     73    GRAPH_TYPEDEFS(typename Graph);
     74
    7575    typedef typename CapacityMap::Value Capacity;
    7676    typedef typename CostMap::Value Cost;
    7777    typedef typename SupplyMap::Value Supply;
    7878
    79     typedef SmartGraph SGraph;
    80     GRAPH_TYPEDEFS(typename SGraph);
    81 
    82     typedef typename SGraph::template EdgeMap<Capacity> SCapacityMap;
    83     typedef typename SGraph::template EdgeMap<Cost> SCostMap;
    84     typedef typename SGraph::template NodeMap<Supply> SSupplyMap;
    85     typedef typename SGraph::template NodeMap<Cost> SPotentialMap;
    86 
    87     typedef typename SGraph::template NodeMap<int> IntNodeMap;
    88     typedef typename SGraph::template NodeMap<bool> BoolNodeMap;
    89     typedef typename SGraph::template NodeMap<Node> NodeNodeMap;
    90     typedef typename SGraph::template NodeMap<Edge> EdgeNodeMap;
    91     typedef typename SGraph::template EdgeMap<int> IntEdgeMap;
    92     typedef typename SGraph::template EdgeMap<bool> BoolEdgeMap;
    93 
    94     typedef typename Graph::template NodeMap<Node> NodeRefMap;
    95     typedef typename Graph::template EdgeMap<Edge> EdgeRefMap;
    96 
    9779    typedef std::vector<Edge> EdgeVector;
     80    typedef std::vector<Node> NodeVector;
     81    typedef std::vector<int> IntVector;
     82    typedef std::vector<bool> BoolVector;
     83    typedef std::vector<Capacity> CapacityVector;
     84    typedef std::vector<Cost> CostVector;
     85    typedef std::vector<Supply> SupplyVector;
     86
     87    typedef typename Graph::template NodeMap<int> IntNodeMap;
    9888
    9989  public:
     
    117107  private:
    118108
    119     /// \brief Map adaptor class for handling reduced edge costs.
    120     ///
    121     /// Map adaptor class for handling reduced edge costs.
    122     class ReducedCostMap : public MapBase<Edge, Cost>
    123     {
    124     private:
    125 
    126       const SGraph &_gr;
    127       const SCostMap &_cost_map;
    128       const SPotentialMap &_pot_map;
    129 
    130     public:
    131 
    132       ///\e
    133       ReducedCostMap( const SGraph &gr,
    134                       const SCostMap &cost_map,
    135                       const SPotentialMap &pot_map ) :
    136         _gr(gr), _cost_map(cost_map), _pot_map(pot_map) {}
    137 
    138       ///\e
    139       Cost operator[](const Edge &e) const {
    140         return _cost_map[e] + _pot_map[_gr.source(e)]
    141                             - _pot_map[_gr.target(e)];
    142       }
    143 
    144     }; //class ReducedCostMap
    145 
    146   private:
    147 
    148109    /// \brief Implementation of the "First Eligible" pivot rule for the
    149110    /// \ref NetworkSimplex "network simplex" algorithm.
     
    158119
    159120      // References to the NetworkSimplex class
    160       NetworkSimplex &_ns;
    161       EdgeVector &_edges;
    162 
     121      const EdgeVector &_edge;
     122      const IntVector  &_source;
     123      const IntVector  &_target;
     124      const CostVector &_cost;
     125      const IntVector  &_state;
     126      const CostVector &_pi;
     127      int &_in_edge;
     128      int _edge_num;
     129
     130      // Pivot rule data
    163131      int _next_edge;
    164132
     
    166134
    167135      /// Constructor
    168       FirstEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
    169         _ns(ns), _edges(edges), _next_edge(0) {}
     136      FirstEligiblePivotRule(NetworkSimplex &ns) :
     137        _edge(ns._edge), _source(ns._source), _target(ns._target),
     138        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
     139        _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
     140      {}
    170141
    171142      /// Find next entering edge
    172143      bool findEnteringEdge() {
    173         Edge e;
    174         for (int i = _next_edge; i < int(_edges.size()); ++i) {
    175           e = _edges[i];
    176           if (_ns._state[e] * _ns._red_cost[e] < 0) {
    177             _ns._in_edge = e;
    178             _next_edge = i + 1;
     144        Cost c;
     145        for (int e = _next_edge; e < _edge_num; ++e) {
     146          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     147          if (c < 0) {
     148            _in_edge = e;
     149            _next_edge = e + 1;
    179150            return true;
    180151          }
    181152        }
    182         for (int i = 0; i < _next_edge; ++i) {
    183           e = _edges[i];
    184           if (_ns._state[e] * _ns._red_cost[e] < 0) {
    185             _ns._in_edge = e;
    186             _next_edge = i + 1;
     153        for (int e = 0; e < _next_edge; ++e) {
     154          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     155          if (c < 0) {
     156            _in_edge = e;
     157            _next_edge = e + 1;
    187158            return true;
    188159          }
     
    190161        return false;
    191162      }
     163
    192164    }; //class FirstEligiblePivotRule
     165
    193166
    194167    /// \brief Implementation of the "Best Eligible" pivot rule for the
     
    204177
    205178      // References to the NetworkSimplex class
    206       NetworkSimplex &_ns;
    207       EdgeVector &_edges;
     179      const EdgeVector &_edge;
     180      const IntVector  &_source;
     181      const IntVector  &_target;
     182      const CostVector &_cost;
     183      const IntVector  &_state;
     184      const CostVector &_pi;
     185      int &_in_edge;
     186      int _edge_num;
    208187
    209188    public:
    210189
    211190      /// Constructor
    212       BestEligiblePivotRule(NetworkSimplex &ns, EdgeVector &edges) :
    213         _ns(ns), _edges(edges) {}
     191      BestEligiblePivotRule(NetworkSimplex &ns) :
     192        _edge(ns._edge), _source(ns._source), _target(ns._target),
     193        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
     194        _in_edge(ns._in_edge), _edge_num(ns._edge_num)
     195      {}
    214196
    215197      /// Find next entering edge
    216198      bool findEnteringEdge() {
    217         Cost min = 0;
    218         Edge e;
    219         for (int i = 0; i < int(_edges.size()); ++i) {
    220           e = _edges[i];
    221           if (_ns._state[e] * _ns._red_cost[e] < min) {
    222             min = _ns._state[e] * _ns._red_cost[e];
    223             _ns._in_edge = e;
     199        Cost c, min = 0;
     200        for (int e = 0; e < _edge_num; ++e) {
     201          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     202          if (c < min) {
     203            min = c;
     204            _in_edge = e;
    224205          }
    225206        }
    226207        return min < 0;
    227208      }
     209
    228210    }; //class BestEligiblePivotRule
     211
    229212
    230213    /// \brief Implementation of the "Block Search" pivot rule for the
     
    240223
    241224      // References to the NetworkSimplex class
    242       NetworkSimplex &_ns;
    243       EdgeVector &_edges;
    244 
     225      const EdgeVector &_edge;
     226      const IntVector  &_source;
     227      const IntVector  &_target;
     228      const CostVector &_cost;
     229      const IntVector  &_state;
     230      const CostVector &_pi;
     231      int &_in_edge;
     232      int _edge_num;
     233
     234      // Pivot rule data
    245235      int _block_size;
    246       int _next_edge, _min_edge;
     236      int _next_edge;
    247237
    248238    public:
    249239
    250240      /// Constructor
    251       BlockSearchPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
    252         _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
     241      BlockSearchPivotRule(NetworkSimplex &ns) :
     242        _edge(ns._edge), _source(ns._source), _target(ns._target),
     243        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
     244        _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
    253245      {
    254246        // The main parameters of the pivot rule
     
    256248        const int MIN_BLOCK_SIZE = 10;
    257249
    258         _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
     250        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
    259251                                MIN_BLOCK_SIZE );
    260252      }
     
    262254      /// Find next entering edge
    263255      bool findEnteringEdge() {
    264         Cost curr, min = 0;
    265         Edge e;
     256        Cost c, min = 0;
    266257        int cnt = _block_size;
    267         int i;
    268         for (i = _next_edge; i < int(_edges.size()); ++i) {
    269           e = _edges[i];
    270           if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
    271             min = curr;
    272             _min_edge = i;
     258        int e, min_edge = _next_edge;
     259        for (e = _next_edge; e < _edge_num; ++e) {
     260          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     261          if (c < min) {
     262            min = c;
     263            min_edge = e;
    273264          }
    274265          if (--cnt == 0) {
     
    278269        }
    279270        if (min == 0 || cnt > 0) {
    280           for (i = 0; i < _next_edge; ++i) {
    281             e = _edges[i];
    282             if ((curr = _ns._state[e] * _ns._red_cost[e]) < min) {
    283               min = curr;
    284               _min_edge = i;
     271          for (e = 0; e < _next_edge; ++e) {
     272            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     273            if (c < min) {
     274              min = c;
     275              min_edge = e;
    285276            }
    286277            if (--cnt == 0) {
     
    291282        }
    292283        if (min >= 0) return false;
    293         _ns._in_edge = _edges[_min_edge];
    294         _next_edge = i;
     284        _in_edge = min_edge;
     285        _next_edge = e;
    295286        return true;
    296287      }
     288
    297289    }; //class BlockSearchPivotRule
     290
    298291
    299292    /// \brief Implementation of the "Candidate List" pivot rule for the
     
    309302
    310303      // References to the NetworkSimplex class
    311       NetworkSimplex &_ns;
    312       EdgeVector &_edges;
    313 
    314       EdgeVector _candidates;
     304      const EdgeVector &_edge;
     305      const IntVector  &_source;
     306      const IntVector  &_target;
     307      const CostVector &_cost;
     308      const IntVector  &_state;
     309      const CostVector &_pi;
     310      int &_in_edge;
     311      int _edge_num;
     312
     313      // Pivot rule data
     314      IntVector _candidates;
    315315      int _list_length, _minor_limit;
    316316      int _curr_length, _minor_count;
    317       int _next_edge, _min_edge;
     317      int _next_edge;
    318318
    319319    public:
    320320
    321321      /// Constructor
    322       CandidateListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
    323         _ns(ns), _edges(edges), _next_edge(0), _min_edge(0)
     322      CandidateListPivotRule(NetworkSimplex &ns) :
     323        _edge(ns._edge), _source(ns._source), _target(ns._target),
     324        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
     325        _in_edge(ns._in_edge), _edge_num(ns._edge_num), _next_edge(0)
    324326      {
    325327        // The main parameters of the pivot rule
     
    329331        const int MIN_MINOR_LIMIT = 3;
    330332
    331         _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edges.size())),
     333        _list_length = std::max( int(LIST_LENGTH_FACTOR * sqrt(_edge_num)),
    332334                                 MIN_LIST_LENGTH );
    333335        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
     
    339341      /// Find next entering edge
    340342      bool findEnteringEdge() {
    341         Cost min, curr;
     343        Cost min, c;
     344        int e, min_edge = _next_edge;
    342345        if (_curr_length > 0 && _minor_count < _minor_limit) {
    343346          // Minor iteration: select the best eligible edge from the
    344347          // current candidate list
    345348          ++_minor_count;
    346           Edge e;
    347349          min = 0;
    348350          for (int i = 0; i < _curr_length; ++i) {
    349351            e = _candidates[i];
    350             curr = _ns._state[e] * _ns._red_cost[e];
    351             if (curr < min) {
    352               min = curr;
    353               _ns._in_edge = e;
     352            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     353            if (c < min) {
     354              min = c;
     355              min_edge = e;
    354356            }
    355             if (curr >= 0) {
     357            if (c >= 0) {
    356358              _candidates[i--] = _candidates[--_curr_length];
    357359            }
    358360          }
    359           if (min < 0) return true;
     361          if (min < 0) {
     362            _in_edge = min_edge;
     363            return true;
     364          }
    360365        }
    361366
    362367        // Major iteration: build a new candidate list
    363         Edge e;
    364368        min = 0;
    365369        _curr_length = 0;
    366         int i;
    367         for (i = _next_edge; i < int(_edges.size()); ++i) {
    368           e = _edges[i];
    369           if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
     370        for (e = _next_edge; e < _edge_num; ++e) {
     371          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     372          if (c < 0) {
    370373            _candidates[_curr_length++] = e;
    371             if (curr < min) {
    372               min = curr;
    373               _min_edge = i;
     374            if (c < min) {
     375              min = c;
     376              min_edge = e;
    374377            }
    375378            if (_curr_length == _list_length) break;
     
    377380        }
    378381        if (_curr_length < _list_length) {
    379           for (i = 0; i < _next_edge; ++i) {
    380             e = _edges[i];
    381             if ((curr = _ns._state[e] * _ns._red_cost[e]) < 0) {
     382          for (e = 0; e < _next_edge; ++e) {
     383            c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     384            if (c < 0) {
    382385              _candidates[_curr_length++] = e;
    383               if (curr < min) {
    384                 min = curr;
    385                 _min_edge = i;
     386              if (c < min) {
     387                min = c;
     388                min_edge = e;
    386389              }
    387390              if (_curr_length == _list_length) break;
     
    391394        if (_curr_length == 0) return false;
    392395        _minor_count = 1;
    393         _ns._in_edge = _edges[_min_edge];
    394         _next_edge = i;
     396        _in_edge = min_edge;
     397        _next_edge = e;
    395398        return true;
    396399      }
     400
    397401    }; //class CandidateListPivotRule
     402
    398403
    399404    /// \brief Implementation of the "Altering Candidate List" pivot rule
     
    409414
    410415      // References to the NetworkSimplex class
    411       NetworkSimplex &_ns;
    412       EdgeVector &_edges;
    413 
    414       EdgeVector _candidates;
    415       SCostMap _cand_cost;
     416      const EdgeVector &_edge;
     417      const IntVector  &_source;
     418      const IntVector  &_target;
     419      const CostVector &_cost;
     420      const IntVector  &_state;
     421      const CostVector &_pi;
     422      int &_in_edge;
     423      int _edge_num;
     424
    416425      int _block_size, _head_length, _curr_length;
    417426      int _next_edge;
     427      IntVector _candidates;
     428      CostVector _cand_cost;
    418429
    419430      // Functor class to compare edges during sort of the candidate list
     
    421432      {
    422433      private:
    423         const SCostMap &_map;
     434        const CostVector &_map;
    424435      public:
    425         SortFunc(const SCostMap &map) : _map(map) {}
    426         bool operator()(const Edge &e1, const Edge &e2) {
    427           return _map[e1] > _map[e2];
     436        SortFunc(const CostVector &map) : _map(map) {}
     437        bool operator()(int left, int right) {
     438          return _map[left] > _map[right];
    428439        }
    429440      };
     
    434445
    435446      /// Constructor
    436       AlteringListPivotRule(NetworkSimplex &ns, EdgeVector &edges) :
    437         _ns(ns), _edges(edges), _cand_cost(_ns._graph),
    438         _next_edge(0), _sort_func(_cand_cost)
     447      AlteringListPivotRule(NetworkSimplex &ns) :
     448        _edge(ns._edge), _source(ns._source), _target(ns._target),
     449        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
     450        _in_edge(ns._in_edge), _edge_num(ns._edge_num),
     451         _next_edge(0), _cand_cost(ns._edge_num),_sort_func(_cand_cost)
    439452      {
    440453        // The main parameters of the pivot rule
     
    444457        const int MIN_HEAD_LENGTH = 3;
    445458
    446         _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edges.size())),
     459        _block_size = std::max( int(BLOCK_SIZE_FACTOR * sqrt(_edge_num)),
    447460                                MIN_BLOCK_SIZE );
    448461        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
     
    455468      bool findEnteringEdge() {
    456469        // Check the current candidate list
    457         Edge e;
    458         for (int ix = 0; ix < _curr_length; ++ix) {
    459           e = _candidates[ix];
    460           if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) >= 0) {
    461             _candidates[ix--] = _candidates[--_curr_length];
     470        int e;
     471        for (int i = 0; i < _curr_length; ++i) {
     472          e = _candidates[i];
     473          _cand_cost[e] = _state[e] *
     474            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     475          if (_cand_cost[e] >= 0) {
     476            _candidates[i--] = _candidates[--_curr_length];
    462477          }
    463478        }
     
    467482        int last_edge = 0;
    468483        int limit = _head_length;
    469         for (int i = _next_edge; i < int(_edges.size()); ++i) {
    470           e = _edges[i];
    471           if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
     484       
     485        for (int e = _next_edge; e < _edge_num; ++e) {
     486          _cand_cost[e] = _state[e] *
     487            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     488          if (_cand_cost[e] < 0) {
    472489            _candidates[_curr_length++] = e;
    473             last_edge = i;
     490            last_edge = e;
    474491          }
    475492          if (--cnt == 0) {
     
    480497        }
    481498        if (_curr_length <= limit) {
    482           for (int i = 0; i < _next_edge; ++i) {
    483             e = _edges[i];
    484             if ((_cand_cost[e] = _ns._state[e] * _ns._red_cost[e]) < 0) {
     499          for (int e = 0; e < _next_edge; ++e) {
     500            _cand_cost[e] = _state[e] *
     501              (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     502            if (_cand_cost[e] < 0) {
    485503              _candidates[_curr_length++] = e;
    486               last_edge = i;
     504              last_edge = e;
    487505            }
    488506            if (--cnt == 0) {
     
    501519
    502520        // Pop the first element of the heap
    503         _ns._in_edge = _candidates[0];
     521        _in_edge = _candidates[0];
    504522        pop_heap( _candidates.begin(), _candidates.begin() + _curr_length,
    505523                  _sort_func );
     
    507525        return true;
    508526      }
     527
    509528    }; //class AlteringListPivotRule
    510529
     
    520539  private:
    521540
    522     // The directed graph the algorithm runs on
    523     SGraph _graph;
    524541    // The original graph
    525     const Graph &_graph_ref;
     542    const Graph &_orig_graph;
    526543    // The original lower bound map
    527     const LowerMap *_lower;
    528     // The capacity map
    529     SCapacityMap _capacity;
    530     // The cost map
    531     SCostMap _cost;
    532     // The supply map
    533     SSupplyMap _supply;
    534     bool _valid_supply;
    535 
    536     // Edge map of the current flow
    537     SCapacityMap _flow;
    538     // Node map of the current potentials
    539     SPotentialMap _potential;
    540 
    541     // The depth node map of the spanning tree structure
    542     IntNodeMap _depth;
    543     // The parent node map of the spanning tree structure
    544     NodeNodeMap _parent;
    545     // The pred_edge node map of the spanning tree structure
    546     EdgeNodeMap _pred_edge;
    547     // The thread node map of the spanning tree structure
    548     NodeNodeMap _thread;
    549     // The forward node map of the spanning tree structure
    550     BoolNodeMap _forward;
    551     // The state edge map
    552     IntEdgeMap _state;
    553     // The artificial root node of the spanning tree structure
    554     Node _root;
    555 
    556     // The reduced cost map
    557     ReducedCostMap _red_cost;
    558 
    559     // The non-artifical edges
    560     EdgeVector _edges;
    561 
    562     // Members for handling the original graph
     544    const LowerMap *_orig_lower;
     545    // The original capacity map
     546    const CapacityMap &_orig_cap;
     547    // The original cost map
     548    const CostMap &_orig_cost;
     549    // The original supply map
     550    const SupplyMap *_orig_supply;
     551    // The source node (if no supply map was given)
     552    Node _orig_source;
     553    // The target node (if no supply map was given)
     554    Node _orig_target;
     555    // The flow value (if no supply map was given)
     556    Capacity _orig_flow_value;
     557
     558    // The flow result map
    563559    FlowMap *_flow_result;
     560    // The potential result map
    564561    PotentialMap *_potential_result;
     562    // Indicate if the flow result map is local
    565563    bool _local_flow;
     564    // Indicate if the potential result map is local
    566565    bool _local_potential;
    567     NodeRefMap _node_ref;
    568     EdgeRefMap _edge_ref;
     566
     567    // The edge references
     568    EdgeVector _edge;
     569    // The node references
     570    NodeVector _node;
     571    // The node ids
     572    IntNodeMap _node_id;
     573    // The source nodes of the edges
     574    IntVector _source;
     575    // The target nodess of the edges
     576    IntVector _target;
     577
     578    // The (modified) capacity vector
     579    CapacityVector _cap;
     580    // The cost vector
     581    CostVector _cost;
     582    // The (modified) supply vector
     583    CostVector _supply;
     584    // The current flow vector
     585    CapacityVector _flow;
     586    // The current potential vector
     587    CostVector _pi;
     588
     589    // The number of nodes in the original graph
     590    int _node_num;
     591    // The number of edges in the original graph
     592    int _edge_num;
     593
     594    // The depth vector of the spanning tree structure
     595    IntVector _depth;
     596    // The parent vector of the spanning tree structure
     597    IntVector _parent;
     598    // The pred_edge vector of the spanning tree structure
     599    IntVector _pred;
     600    // The thread vector of the spanning tree structure
     601    IntVector _thread;
     602    // The forward vector of the spanning tree structure
     603    BoolVector _forward;
     604    // The state vector
     605    IntVector _state;
     606    // The root node
     607    int _root;
    569608
    570609    // The entering edge of the current pivot iteration
    571     Edge _in_edge;
     610    int _in_edge;
    572611
    573612    // Temporary nodes used in the current pivot iteration
    574     Node join, u_in, v_in, u_out, v_out;
    575     Node right, first, second, last;
    576     Node stem, par_stem, new_stem;
     613    int join, u_in, v_in, u_out, v_out;
     614    int right, first, second, last;
     615    int stem, par_stem, new_stem;
     616
    577617    // The maximum augment amount along the found cycle in the current
    578618    // pivot iteration
    579619    Capacity delta;
    580620
    581   public :
     621  public:
    582622
    583623    /// \brief General constructor (with lower bounds).
     
    595635                    const CostMap &cost,
    596636                    const SupplyMap &supply ) :
    597       _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
    598       _cost(_graph), _supply(_graph), _flow(_graph),
    599       _potential(_graph), _depth(_graph), _parent(_graph),
    600       _pred_edge(_graph), _thread(_graph), _forward(_graph),
    601       _state(_graph), _red_cost(_graph, _cost, _potential),
     637      _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
     638      _orig_cost(cost), _orig_supply(&supply),
    602639      _flow_result(NULL), _potential_result(NULL),
    603640      _local_flow(false), _local_potential(false),
    604       _node_ref(graph), _edge_ref(graph)
    605     {
    606       // Check the sum of supply values
    607       Supply sum = 0;
    608       for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
    609         sum += supply[n];
    610       if (!(_valid_supply = sum == 0)) return;
    611 
    612       // Copy _graph_ref to _graph
    613       _graph.reserveNode(countNodes(_graph_ref) + 1);
    614       _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
    615       copyGraph(_graph, _graph_ref)
    616         .edgeMap(_capacity, capacity)
    617         .edgeMap(_cost, cost)
    618         .nodeMap(_supply, supply)
    619         .nodeRef(_node_ref)
    620         .edgeRef(_edge_ref)
    621         .run();
    622 
    623       // Remove non-zero lower bounds
    624       for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
    625         if (lower[e] != 0) {
    626         _capacity[_edge_ref[e]] = capacity[e] - lower[e];
    627           _supply[_node_ref[_graph_ref.source(e)]] -= lower[e];
    628           _supply[_node_ref[_graph_ref.target(e)]] += lower[e];
    629       }
    630       }
    631     }
     641      _node_id(graph)
     642    {}
    632643
    633644    /// \brief General constructor (without lower bounds).
     
    643654                    const CostMap &cost,
    644655                    const SupplyMap &supply ) :
    645       _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
    646       _cost(_graph), _supply(_graph), _flow(_graph),
    647       _potential(_graph), _depth(_graph), _parent(_graph),
    648       _pred_edge(_graph), _thread(_graph), _forward(_graph),
    649       _state(_graph), _red_cost(_graph, _cost, _potential),
     656      _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
     657      _orig_cost(cost), _orig_supply(&supply),
    650658      _flow_result(NULL), _potential_result(NULL),
    651659      _local_flow(false), _local_potential(false),
    652       _node_ref(graph), _edge_ref(graph)
    653     {
    654       // Check the sum of supply values
    655       Supply sum = 0;
    656       for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
    657         sum += supply[n];
    658       if (!(_valid_supply = sum == 0)) return;
    659 
    660       // Copy _graph_ref to _graph
    661       _graph.reserveNode(countNodes(_graph_ref) + 1);
    662       _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
    663       copyGraph(_graph, _graph_ref)
    664         .edgeMap(_capacity, capacity)
    665         .edgeMap(_cost, cost)
    666         .nodeMap(_supply, supply)
    667         .nodeRef(_node_ref)
    668         .edgeRef(_edge_ref)
    669         .run();
    670     }
     660      _node_id(graph)
     661    {}
    671662
    672663    /// \brief Simple constructor (with lower bounds).
     
    686677                    const CapacityMap &capacity,
    687678                    const CostMap &cost,
    688                     typename Graph::Node s,
    689                     typename Graph::Node t,
    690                     typename SupplyMap::Value flow_value ) :
    691       _graph(), _graph_ref(graph), _lower(&lower), _capacity(_graph),
    692       _cost(_graph), _supply(_graph, 0), _flow(_graph),
    693       _potential(_graph), _depth(_graph), _parent(_graph),
    694       _pred_edge(_graph), _thread(_graph), _forward(_graph),
    695       _state(_graph), _red_cost(_graph, _cost, _potential),
     679                    Node s, Node t,
     680                    Capacity flow_value ) :
     681      _orig_graph(graph), _orig_lower(&lower), _orig_cap(capacity),
     682      _orig_cost(cost), _orig_supply(NULL),
     683      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
    696684      _flow_result(NULL), _potential_result(NULL),
    697685      _local_flow(false), _local_potential(false),
    698       _node_ref(graph), _edge_ref(graph)
    699     {
    700       // Copy _graph_ref to graph
    701       _graph.reserveNode(countNodes(_graph_ref) + 1);
    702       _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
    703       copyGraph(_graph, _graph_ref)
    704         .edgeMap(_capacity, capacity)
    705         .edgeMap(_cost, cost)
    706         .nodeRef(_node_ref)
    707         .edgeRef(_edge_ref)
    708         .run();
    709 
    710       // Remove non-zero lower bounds
    711       _supply[_node_ref[s]] =  flow_value;
    712       _supply[_node_ref[t]] = -flow_value;
    713       for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e) {
    714         if (lower[e] != 0) {
    715         _capacity[_edge_ref[e]] = capacity[e] - lower[e];
    716           _supply[_node_ref[_graph_ref.source(e)]] -= lower[e];
    717           _supply[_node_ref[_graph_ref.target(e)]] += lower[e];
    718       }
    719       }
    720       _valid_supply = true;
    721     }
     686      _node_id(graph)
     687    {}
    722688
    723689    /// \brief Simple constructor (without lower bounds).
     
    735701                    const CapacityMap &capacity,
    736702                    const CostMap &cost,
    737                     typename Graph::Node s,
    738                     typename Graph::Node t,
    739                     typename SupplyMap::Value flow_value ) :
    740       _graph(), _graph_ref(graph), _lower(NULL), _capacity(_graph),
    741       _cost(_graph), _supply(_graph, 0), _flow(_graph),
    742       _potential(_graph), _depth(_graph), _parent(_graph),
    743       _pred_edge(_graph), _thread(_graph), _forward(_graph),
    744       _state(_graph), _red_cost(_graph, _cost, _potential),
     703                    Node s, Node t,
     704                    Capacity flow_value ) :
     705      _orig_graph(graph), _orig_lower(NULL), _orig_cap(capacity),
     706      _orig_cost(cost), _orig_supply(NULL),
     707      _orig_source(s), _orig_target(t), _orig_flow_value(flow_value),
    745708      _flow_result(NULL), _potential_result(NULL),
    746709      _local_flow(false), _local_potential(false),
    747       _node_ref(graph), _edge_ref(graph)
    748     {
    749       // Copy _graph_ref to graph
    750       _graph.reserveNode(countNodes(_graph_ref) + 1);
    751       _graph.reserveEdge(countEdges(_graph_ref) + countNodes(_graph_ref));
    752       copyGraph(_graph, _graph_ref)
    753         .edgeMap(_capacity, capacity)
    754         .edgeMap(_cost, cost)
    755         .nodeRef(_node_ref)
    756         .edgeRef(_edge_ref)
    757         .run();
    758       _supply[_node_ref[s]] =  flow_value;
    759       _supply[_node_ref[t]] = -flow_value;
    760       _valid_supply = true;
    761     }
     710      _node_id(graph)
     711    {}
    762712
    763713    /// Destructor.
     
    895845    Cost totalCost() const {
    896846      Cost c = 0;
    897       for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
    898         c += (*_flow_result)[e] * _cost[_edge_ref[e]];
     847      for (EdgeIt e(_orig_graph); e != INVALID; ++e)
     848        c += (*_flow_result)[e] * _orig_cost[e];
    899849      return c;
    900850    }
     
    904854  private:
    905855
    906     // Extend the underlying graph and initialize all the node and edge maps
     856    // Initialize internal data structures
    907857    bool init() {
    908       if (!_valid_supply) return false;
    909 
    910858      // Initialize result maps
    911859      if (!_flow_result) {
    912         _flow_result = new FlowMap(_graph_ref);
     860        _flow_result = new FlowMap(_orig_graph);
    913861        _local_flow = true;
    914862      }
    915863      if (!_potential_result) {
    916         _potential_result = new PotentialMap(_graph_ref);
     864        _potential_result = new PotentialMap(_orig_graph);
    917865        _local_potential = true;
    918866      }
    919 
    920       // Initialize the edge vector and the edge maps
    921       _edges.reserve(countEdges(_graph));
    922       for (EdgeIt e(_graph); e != INVALID; ++e) {
    923         _edges.push_back(e);
    924         _flow[e] = 0;
    925         _state[e] = STATE_LOWER;
    926       }
    927 
    928       // Add an artificial root node to the graph
    929       _root = _graph.addNode();
    930       _parent[_root] = INVALID;
    931       _pred_edge[_root] = INVALID;
     867     
     868      // Initialize vectors
     869      _node_num = countNodes(_orig_graph);
     870      _edge_num = countEdges(_orig_graph);
     871      int all_node_num = _node_num + 1;
     872      int all_edge_num = _edge_num + _node_num;
     873     
     874      _edge.resize(_edge_num);
     875      _node.reserve(_node_num);
     876      _source.resize(all_edge_num);
     877      _target.resize(all_edge_num);
     878     
     879      _cap.resize(all_edge_num);
     880      _cost.resize(all_edge_num);
     881      _supply.resize(all_node_num);
     882      _flow.resize(all_edge_num, 0);
     883      _pi.resize(all_node_num, 0);
     884     
     885      _depth.resize(all_node_num);
     886      _parent.resize(all_node_num);
     887      _pred.resize(all_node_num);
     888      _thread.resize(all_node_num);
     889      _forward.resize(all_node_num);
     890      _state.resize(all_edge_num, STATE_LOWER);
     891     
     892      // Initialize node related data
     893      bool valid_supply = true;
     894      if (_orig_supply) {
     895        Supply sum = 0;
     896        int i = 0;
     897        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
     898          _node.push_back(n);
     899          _node_id[n] = i;
     900          _supply[i] = (*_orig_supply)[n];
     901          sum += _supply[i];
     902        }
     903        valid_supply = (sum == 0);
     904      } else {
     905        int i = 0;
     906        for (NodeIt n(_orig_graph); n != INVALID; ++n, ++i) {
     907          _node.push_back(n);
     908          _node_id[n] = i;
     909          _supply[i] = 0;
     910        }
     911        _supply[_node_id[_orig_source]] =  _orig_flow_value;
     912        _supply[_node_id[_orig_target]] = -_orig_flow_value;
     913      }
     914      if (!valid_supply) return false;
     915     
     916      // Set data for the artificial root node
     917      _root = _node_num;
    932918      _depth[_root] = 0;
     919      _parent[_root] = -1;
     920      _pred[_root] = -1;
     921      _thread[_root] = 0;
    933922      _supply[_root] = 0;
    934       _potential[_root] = 0;
    935 
    936       // Add artificial edges to the graph and initialize the node maps
    937       // of the spanning tree data structure
    938       Node last = _root;
    939       Edge e;
     923      _pi[_root] = 0;
     924     
     925      // Store the edges in a mixed order
     926      int k = std::max(int(sqrt(_edge_num)), 10);
     927      int i = 0;
     928      for (EdgeIt e(_orig_graph); e != INVALID; ++e) {
     929        _edge[i] = e;
     930        if ((i += k) >= _edge_num) i = (i % k) + 1;
     931      }
     932
     933      // Initialize edge maps
     934      for (int i = 0; i != _edge_num; ++i) {
     935        Edge e = _edge[i];
     936        _source[i] = _node_id[_orig_graph.source(e)];
     937        _target[i] = _node_id[_orig_graph.target(e)];
     938        _cost[i] = _orig_cost[e];
     939        _cap[i] = _orig_cap[e];
     940      }
     941
     942      // Remove non-zero lower bounds
     943      if (_orig_lower) {
     944        for (int i = 0; i != _edge_num; ++i) {
     945          Capacity c = (*_orig_lower)[_edge[i]];
     946          if (c != 0) {
     947            _cap[i] -= c;
     948            _supply[_source[i]] -= c;
     949            _supply[_target[i]] += c;
     950          }
     951        }
     952      }
     953
     954      // Add artificial edges and initialize the spanning tree data structure
    940955      Cost max_cost = std::numeric_limits<Cost>::max() / 4;
    941       for (NodeIt u(_graph); u != INVALID; ++u) {
    942         if (u == _root) continue;
    943         _thread[last] = u;
    944         last = u;
     956      Capacity max_cap = std::numeric_limits<Capacity>::max();
     957      for (int u = 0, e = _edge_num; u != _node_num; ++u, ++e) {
     958        _thread[u] = u + 1;
     959        _depth[u] = 1;
    945960        _parent[u] = _root;
    946         _depth[u] = 1;
     961        _pred[u] = e;
    947962        if (_supply[u] >= 0) {
    948           e = _graph.addEdge(u, _root);
    949963          _flow[e] = _supply[u];
    950964          _forward[u] = true;
    951           _potential[u] = -max_cost;
     965          _pi[u] = -max_cost;
    952966        } else {
    953           e = _graph.addEdge(_root, u);
    954967          _flow[e] = -_supply[u];
    955968          _forward[u] = false;
    956           _potential[u] = max_cost;
     969          _pi[u] = max_cost;
    957970        }
    958971        _cost[e] = max_cost;
    959         _capacity[e] = std::numeric_limits<Capacity>::max();
     972        _cap[e] = max_cap;
    960973        _state[e] = STATE_TREE;
    961         _pred_edge[u] = e;
    962       }
    963       _thread[last] = _root;
     974      }
    964975
    965976      return true;
     
    968979    // Find the join node
    969980    void findJoinNode() {
    970       Node u = _graph.source(_in_edge);
    971       Node v = _graph.target(_in_edge);
     981      int u = _source[_in_edge];
     982      int v = _target[_in_edge];
     983      while (_depth[u] > _depth[v]) u = _parent[u];
     984      while (_depth[v] > _depth[u]) v = _parent[v];
    972985      while (u != v) {
    973         if (_depth[u] == _depth[v]) {
    974           u = _parent[u];
    975           v = _parent[v];
    976         }
    977         else if (_depth[u] > _depth[v]) u = _parent[u];
    978         else v = _parent[v];
     986        u = _parent[u];
     987        v = _parent[v];
    979988      }
    980989      join = u;
    981990    }
    982991
    983     // Find the leaving edge of the cycle and returns true if the 
     992    // Find the leaving edge of the cycle and returns true if the
    984993    // leaving edge is not the same as the entering edge
    985994    bool findLeavingEdge() {
     
    987996      // of the cycle
    988997      if (_state[_in_edge] == STATE_LOWER) {
    989         first  = _graph.source(_in_edge);
    990         second = _graph.target(_in_edge);
     998        first  = _source[_in_edge];
     999        second = _target[_in_edge];
    9911000      } else {
    992         first  = _graph.target(_in_edge);
    993         second = _graph.source(_in_edge);
    994       }
    995       delta = _capacity[_in_edge];
    996       bool result = false;
     1001        first  = _target[_in_edge];
     1002        second = _source[_in_edge];
     1003      }
     1004      delta = _cap[_in_edge];
     1005      int result = 0;
    9971006      Capacity d;
    998       Edge e;
     1007      int e;
    9991008
    10001009      // Search the cycle along the path form the first node to the root
    1001       for (Node u = first; u != join; u = _parent[u]) {
    1002         e = _pred_edge[u];
    1003         d = _forward[u] ? _flow[e] : _capacity[e] - _flow[e];
     1010      for (int u = first; u != join; u = _parent[u]) {
     1011        e = _pred[u];
     1012        d = _forward[u] ? _flow[e] : _cap[e] - _flow[e];
    10041013        if (d < delta) {
    10051014          delta = d;
    10061015          u_out = u;
    1007           u_in = first;
    1008           v_in = second;
    1009           result = true;
     1016          result = 1;
    10101017        }
    10111018      }
    10121019      // Search the cycle along the path form the second node to the root
    1013       for (Node u = second; u != join; u = _parent[u]) {
    1014         e = _pred_edge[u];
    1015         d = _forward[u] ? _capacity[e] - _flow[e] : _flow[e];
     1020      for (int u = second; u != join; u = _parent[u]) {
     1021        e = _pred[u];
     1022        d = _forward[u] ? _cap[e] - _flow[e] : _flow[e];
    10161023        if (d <= delta) {
    10171024          delta = d;
    10181025          u_out = u;
    1019           u_in = second;
    1020           v_in = first;
    1021           result = true;
    1022         }
    1023       }
    1024       return result;
    1025     }
    1026 
    1027     // Change _flow and _state edge maps
    1028     void changeFlows(bool change) {
     1026          result = 2;
     1027        }
     1028      }
     1029
     1030      if (result == 1) {
     1031        u_in = first;
     1032        v_in = second;
     1033      } else {
     1034        u_in = second;
     1035        v_in = first;
     1036      }
     1037      return result != 0;
     1038    }
     1039
     1040    // Change _flow and _state vectors
     1041    void changeFlow(bool change) {
    10291042      // Augment along the cycle
    10301043      if (delta > 0) {
    10311044        Capacity val = _state[_in_edge] * delta;
    10321045        _flow[_in_edge] += val;
    1033         for (Node u = _graph.source(_in_edge); u != join; u = _parent[u]) {
    1034           _flow[_pred_edge[u]] += _forward[u] ? -val : val;
    1035         }
    1036         for (Node u = _graph.target(_in_edge); u != join; u = _parent[u]) {
    1037           _flow[_pred_edge[u]] += _forward[u] ? val : -val;
     1046        for (int u = _source[_in_edge]; u != join; u = _parent[u]) {
     1047          _flow[_pred[u]] += _forward[u] ? -val : val;
     1048        }
     1049        for (int u = _target[_in_edge]; u != join; u = _parent[u]) {
     1050          _flow[_pred[u]] += _forward[u] ? val : -val;
    10381051        }
    10391052      }
     
    10411054      if (change) {
    10421055        _state[_in_edge] = STATE_TREE;
    1043         _state[_pred_edge[u_out]] =
    1044           (_flow[_pred_edge[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
     1056        _state[_pred[u_out]] =
     1057          (_flow[_pred[u_out]] == 0) ? STATE_LOWER : STATE_UPPER;
    10451058      } else {
    10461059        _state[_in_edge] = -_state[_in_edge];
     
    10481061    }
    10491062
    1050     // Update _thread and _parent node maps
     1063    // Update _thread and _parent vectors
    10511064    void updateThreadParent() {
    1052       Node u;
     1065      int u;
    10531066      v_out = _parent[u_out];
    10541067
     
    11061119    }
    11071120
    1108     // Update _pred_edge and _forward node maps.
     1121    // Update _pred and _forward vectors
    11091122    void updatePredEdge() {
    1110       Node u = u_out, v;
     1123      int u = u_out, v;
    11111124      while (u != u_in) {
    11121125        v = _parent[u];
    1113         _pred_edge[u] = _pred_edge[v];
     1126        _pred[u] = _pred[v];
    11141127        _forward[u] = !_forward[v];
    11151128        u = v;
    11161129      }
    1117       _pred_edge[u_in] = _in_edge;
    1118       _forward[u_in] = (u_in == _graph.source(_in_edge));
    1119     }
    1120 
    1121     // Update _depth and _potential node maps
     1130      _pred[u_in] = _in_edge;
     1131      _forward[u_in] = (u_in == _source[_in_edge]);
     1132    }
     1133
     1134    // Update _depth and _potential vectors
    11221135    void updateDepthPotential() {
    11231136      _depth[u_in] = _depth[v_in] + 1;
    11241137      Cost sigma = _forward[u_in] ?
    1125         _potential[v_in] - _potential[u_in] - _cost[_pred_edge[u_in]] :
    1126         _potential[v_in] - _potential[u_in] + _cost[_pred_edge[u_in]];
    1127       _potential[u_in] += sigma;
    1128       for(Node u = _thread[u_in]; _parent[u] != INVALID; u = _thread[u]) {
     1138        _pi[v_in] - _pi[u_in] - _cost[_pred[u_in]] :
     1139        _pi[v_in] - _pi[u_in] + _cost[_pred[u_in]];
     1140      _pi[u_in] += sigma;
     1141      for(int u = _thread[u_in]; _parent[u] != -1; u = _thread[u]) {
    11291142        _depth[u] = _depth[_parent[u]] + 1;
    11301143        if (_depth[u] <= _depth[u_in]) break;
    1131         _potential[u] += sigma;
     1144        _pi[u] += sigma;
    11321145      }
    11331146    }
     
    11531166    template<class PivotRuleImplementation>
    11541167    bool start() {
    1155       PivotRuleImplementation pivot(*this, _edges);
     1168      PivotRuleImplementation pivot(*this);
    11561169
    11571170      // Execute the network simplex algorithm
     
    11591172        findJoinNode();
    11601173        bool change = findLeavingEdge();
    1161         changeFlows(change);
     1174        changeFlow(change);
    11621175        if (change) {
    11631176          updateThreadParent();
     
    11681181
    11691182      // Check if the flow amount equals zero on all the artificial edges
    1170       for (InEdgeIt e(_graph, _root); e != INVALID; ++e)
     1183      for (int e = _edge_num; e != _edge_num + _node_num; ++e) {
    11711184        if (_flow[e] > 0) return false;
    1172       for (OutEdgeIt e(_graph, _root); e != INVALID; ++e)
    1173         if (_flow[e] > 0) return false;
     1185      }
    11741186
    11751187      // Copy flow values to _flow_result
    1176       if (_lower) {
    1177         for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
    1178           (*_flow_result)[e] = (*_lower)[e] + _flow[_edge_ref[e]];
     1188      if (_orig_lower) {
     1189        for (int i = 0; i != _edge_num; ++i) {
     1190          Edge e = _edge[i];
     1191          (*_flow_result)[e] = (*_orig_lower)[e] + _flow[i];
     1192        }
    11791193      } else {
    1180         for (typename Graph::EdgeIt e(_graph_ref); e != INVALID; ++e)
    1181           (*_flow_result)[e] = _flow[_edge_ref[e]];
     1194        for (int i = 0; i != _edge_num; ++i) {
     1195          (*_flow_result)[_edge[i]] = _flow[i];
     1196        }
    11821197      }
    11831198      // Copy potential values to _potential_result
    1184       for (typename Graph::NodeIt n(_graph_ref); n != INVALID; ++n)
    1185         (*_potential_result)[n] = _potential[_node_ref[n]];
     1199      for (int i = 0; i != _node_num; ++i) {
     1200        (*_potential_result)[_node[i]] = _pi[i];
     1201      }
    11861202
    11871203      return true;
Note: See TracChangeset for help on using the changeset viewer.