COIN-OR::LEMON - Graph Library

Changeset 2625:c51b320bc51c in lemon-0.x for lemon/cost_scaling.h


Ignore:
Timestamp:
10/15/08 14:04:11 (11 years ago)
Author:
Peter Kovacs
Branch:
default
Phase:
public
Convert:
svn:c9d7d8f5-90d6-0310-b91f-818b3a526b0e/lemon/trunk@3510
Message:

Major improvement in the cost scaling algorithm

  • Add a new variant that use the partial augment-relabel method.
  • Use this method instead of push-relabel by default.
  • Use the "Early Termination" heuristic instead of "Price Refinement".

Using the new method and heuristic the algorithm proved to be
2-2.5 times faster on all input files.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/cost_scaling.h

    r2623 r2625  
    2121
    2222/// \ingroup min_cost_flow
    23 ///
    2423/// \file
    2524/// \brief Cost scaling algorithm for finding a minimum cost flow.
     
    4342  ///
    4443  /// \ref CostScaling implements the cost scaling algorithm performing
    45   /// generalized push-relabel operations for finding a minimum cost
     44  /// augment/push and relabel operations for finding a minimum cost
    4645  /// flow.
    4746  ///
     
    117116
    118117      ///\e
    119       typename Map::Value operator[](const ResEdge &e) const {
    120         return ResGraph::forward(e) ?  _cost_map[e] : -_cost_map[e];
     118      inline typename Map::Value operator[](const ResEdge &e) const {
     119        return ResGraph::forward(e) ? _cost_map[e] : -_cost_map[e];
    121120      }
    122121
     
    143142
    144143      ///\e
    145       LCost operator[](const Edge &e) const {
     144      inline LCost operator[](const Edge &e) const {
    146145        return _cost_map[e] + _pot_map[_gr.source(e)]
    147146                            - _pot_map[_gr.target(e)];
     
    149148
    150149    }; //class ReducedCostMap
    151 
    152   private:
    153 
    154     // Scaling factor
    155     static const int ALPHA = 4;
    156 
    157     // Paramters for heuristics
    158     static const int BF_HEURISTIC_EPSILON_BOUND = 5000;
    159     static const int BF_HEURISTIC_BOUND_FACTOR  = 3;
    160150
    161151  private:
     
    192182    // The epsilon parameter used for cost scaling
    193183    LCost _epsilon;
     184    // The scaling factor
     185    int _alpha;
    194186
    195187  public:
     
    214206      _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
    215207    {
    216       // Removing non-zero lower bounds
     208      // Remove non-zero lower bounds
    217209      _capacity = subMap(capacity, lower);
    218210      Supply sum = 0;
     
    246238      _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
    247239    {
    248       // Checking the sum of supply values
     240      // Check the sum of supply values
    249241      Supply sum = 0;
    250242      for (NodeIt n(_graph); n != INVALID; ++n) sum += _supply[n];
     
    275267      _res_graph(NULL), _red_cost(NULL), _excess(graph, 0)
    276268    {
    277       // Removing nonzero lower bounds
     269      // Remove nonzero lower bounds
    278270      _capacity = subMap(capacity, lower);
    279271      for (NodeIt n(_graph); n != INVALID; ++n) {
     
    360352    /// Run the algorithm.
    361353    ///
     354    /// \param partial_augment By default the algorithm performs
     355    /// partial augment and relabel operations in the cost scaling
     356    /// phases. Set this parameter to \c false for using local push and
     357    /// relabel operations instead.
     358    ///
    362359    /// \return \c true if a feasible flow can be found.
    363     bool run() {
    364       return init() && start();
     360    bool run(bool partial_augment = true) {
     361      if (partial_augment) {
     362        return init() && startPartialAugment();
     363      } else {
     364        return init() && startPushRelabel();
     365      }
    365366    }
    366367
     
    434435    bool init() {
    435436      if (!_valid_supply) return false;
    436 
    437       // Initializing flow and potential maps
     437      // The scaling factor
     438      _alpha = 8;
     439
     440      // Initialize flow and potential maps
    438441      if (!_flow) {
    439442        _flow = new FlowMap(_graph);
     
    448451      _res_graph = new ResGraph(_graph, _capacity, *_flow);
    449452
    450       // Initializing the scaled cost map and the epsilon parameter
     453      // Initialize the scaled cost map and the epsilon parameter
    451454      Cost max_cost = 0;
    452455      int node_num = countNodes(_graph);
    453456      for (EdgeIt e(_graph); e != INVALID; ++e) {
    454         _cost[e] = LCost(_orig_cost[e]) * node_num * ALPHA;
     457        _cost[e] = LCost(_orig_cost[e]) * node_num * _alpha;
    455458        if (_orig_cost[e] > max_cost) max_cost = _orig_cost[e];
    456459      }
    457460      _epsilon = max_cost * node_num;
    458461
    459       // Finding a feasible flow using Circulation
     462      // Find a feasible flow using Circulation
    460463      Circulation< Graph, ConstMap<Edge, Capacity>, CapacityEdgeMap,
    461464                   SupplyMap >
     
    465468    }
    466469
    467 
    468     /// Execute the algorithm.
    469     bool start() {
     470    /// Execute the algorithm performing partial augmentation and
     471    /// relabel operations.
     472    bool startPartialAugment() {
     473      // Paramters for heuristics
     474      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
     475      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
     476      // Maximum augment path length
     477      const int MAX_PATH_LENGTH = 4;
     478
     479      // Variables
     480      typename Graph::template NodeMap<Edge> pred_edge(_graph);
     481      typename Graph::template NodeMap<bool> forward(_graph);
     482      typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
     483      typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
     484      typename Graph::template NodeMap<bool> next_dir(_graph);
    470485      std::deque<Node> active_nodes;
    471       typename Graph::template NodeMap<bool> hyper(_graph, false);
     486      std::vector<Node> path_nodes;
    472487
    473488      int node_num = countNodes(_graph);
    474       for ( ; _epsilon >= 1; _epsilon = _epsilon < ALPHA && _epsilon > 1 ?
    475                                         1 : _epsilon / ALPHA )
     489      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
     490                                        1 : _epsilon / _alpha )
    476491      {
    477         // Performing price refinement heuristic using Bellman-Ford
    478         // algorithm
     492        // "Early Termination" heuristic: use Bellman-Ford algorithm
     493        // to check if the current flow is optimal
    479494        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
    480495          typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
    481           ShiftCostMap shift_cost(_res_cost, _epsilon);
     496          ShiftCostMap shift_cost(_res_cost, 1);
    482497          BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
    483498          bf.init(0);
     
    486501          for (int i = 0; i < K && !done; ++i)
    487502            done = bf.processNextWeakRound();
    488           if (done) {
    489             for (NodeIt n(_graph); n != INVALID; ++n)
    490               (*_potential)[n] = bf.dist(n);
    491             continue;
    492           }
    493         }
    494 
    495         // Saturating edges not satisfying the optimality condition
     503          if (done) break;
     504        }
     505
     506        // Saturate edges not satisfying the optimality condition
    496507        Capacity delta;
    497508        for (EdgeIt e(_graph); e != INVALID; ++e) {
     
    509520        }
    510521
    511         // Finding active nodes (i.e. nodes with positive excess)
    512         for (NodeIt n(_graph); n != INVALID; ++n)
     522        // Find active nodes (i.e. nodes with positive excess)
     523        for (NodeIt n(_graph); n != INVALID; ++n) {
    513524          if (_excess[n] > 0) active_nodes.push_back(n);
    514 
    515         // Performing push and relabel operations
     525        }
     526
     527        // Initialize the next edge maps
     528        for (NodeIt n(_graph); n != INVALID; ++n) {
     529          next_out[n] = OutEdgeIt(_graph, n);
     530          next_in[n] = InEdgeIt(_graph, n);
     531          next_dir[n] = true;
     532        }
     533
     534        // Perform partial augment and relabel operations
    516535        while (active_nodes.size() > 0) {
     536          // Select an active node (FIFO selection)
     537          if (_excess[active_nodes[0]] <= 0) {
     538            active_nodes.pop_front();
     539            continue;
     540          }
     541          Node start = active_nodes[0];
     542          path_nodes.clear();
     543          path_nodes.push_back(start);
     544
     545          // Find an augmenting path from the start node
     546          Node u, tip = start;
     547          LCost min_red_cost;
     548          while ( _excess[tip] >= 0 &&
     549                  int(path_nodes.size()) <= MAX_PATH_LENGTH )
     550          {
     551            if (next_dir[tip]) {
     552              for (OutEdgeIt e = next_out[tip]; e != INVALID; ++e) {
     553                if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
     554                  u = _graph.target(e);
     555                  pred_edge[u] = e;
     556                  forward[u] = true;
     557                  next_out[tip] = e;
     558                  tip = u;
     559                  path_nodes.push_back(tip);
     560                  goto next_step;
     561                }
     562              }
     563              next_dir[tip] = false;
     564            }
     565            for (InEdgeIt e = next_in[tip]; e != INVALID; ++e) {
     566              if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
     567                u = _graph.source(e);
     568                pred_edge[u] = e;
     569                forward[u] = false;
     570                next_in[tip] = e;
     571                tip = u;
     572                path_nodes.push_back(tip);
     573                goto next_step;
     574              }
     575            }
     576
     577            // Relabel tip node
     578            min_red_cost = std::numeric_limits<LCost>::max() / 2;
     579            for (OutEdgeIt oe(_graph, tip); oe != INVALID; ++oe) {
     580              if ( _capacity[oe] - (*_flow)[oe] > 0 &&
     581                   (*_red_cost)[oe] < min_red_cost )
     582                min_red_cost = (*_red_cost)[oe];
     583            }
     584            for (InEdgeIt ie(_graph, tip); ie != INVALID; ++ie) {
     585              if ((*_flow)[ie] > 0 && -(*_red_cost)[ie] < min_red_cost)
     586                min_red_cost = -(*_red_cost)[ie];
     587            }
     588            (*_potential)[tip] -= min_red_cost + _epsilon;
     589
     590            // Reset the next edge maps
     591            next_out[tip] = OutEdgeIt(_graph, tip);
     592            next_in[tip] = InEdgeIt(_graph, tip);
     593            next_dir[tip] = true;
     594
     595            // Step back
     596            if (tip != start) {
     597              path_nodes.pop_back();
     598              tip = path_nodes[path_nodes.size()-1];
     599            }
     600
     601          next_step:
     602            continue;
     603          }
     604
     605          // Augment along the found path (as much flow as possible)
     606          Capacity delta;
     607          for (int i = 1; i < int(path_nodes.size()); ++i) {
     608            u = path_nodes[i];
     609            delta = forward[u] ?
     610              _capacity[pred_edge[u]] - (*_flow)[pred_edge[u]] :
     611              (*_flow)[pred_edge[u]];
     612            delta = std::min(delta, _excess[path_nodes[i-1]]);
     613            (*_flow)[pred_edge[u]] += forward[u] ? delta : -delta;
     614            _excess[path_nodes[i-1]] -= delta;
     615            _excess[u] += delta;
     616            if (_excess[u] > 0 && _excess[u] <= delta) active_nodes.push_back(u);
     617          }
     618        }
     619      }
     620
     621      // Compute node potentials for the original costs
     622      ResidualCostMap<CostMap> res_cost(_orig_cost);
     623      BellmanFord< ResGraph, ResidualCostMap<CostMap> >
     624        bf(*_res_graph, res_cost);
     625      bf.init(0); bf.start();
     626      for (NodeIt n(_graph); n != INVALID; ++n)
     627        (*_potential)[n] = bf.dist(n);
     628
     629      // Handle non-zero lower bounds
     630      if (_lower) {
     631        for (EdgeIt e(_graph); e != INVALID; ++e)
     632          (*_flow)[e] += (*_lower)[e];
     633      }
     634      return true;
     635    }
     636
     637    /// Execute the algorithm performing push and relabel operations.
     638    bool startPushRelabel() {
     639      // Paramters for heuristics
     640      const int BF_HEURISTIC_EPSILON_BOUND = 1000;
     641      const int BF_HEURISTIC_BOUND_FACTOR  = 3;
     642
     643      typename Graph::template NodeMap<bool> hyper(_graph, false);
     644      typename Graph::template NodeMap<Edge> pred_edge(_graph);
     645      typename Graph::template NodeMap<bool> forward(_graph);
     646      typename Graph::template NodeMap<OutEdgeIt> next_out(_graph);
     647      typename Graph::template NodeMap<InEdgeIt> next_in(_graph);
     648      typename Graph::template NodeMap<bool> next_dir(_graph);
     649      std::deque<Node> active_nodes;
     650
     651      int node_num = countNodes(_graph);
     652      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
     653                                        1 : _epsilon / _alpha )
     654      {
     655        // "Early Termination" heuristic: use Bellman-Ford algorithm
     656        // to check if the current flow is optimal
     657        if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
     658          typedef ShiftMap< ResidualCostMap<LargeCostMap> > ShiftCostMap;
     659          ShiftCostMap shift_cost(_res_cost, 1);
     660          BellmanFord<ResGraph, ShiftCostMap> bf(*_res_graph, shift_cost);
     661          bf.init(0);
     662          bool done = false;
     663          int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(node_num));
     664          for (int i = 0; i < K && !done; ++i)
     665            done = bf.processNextWeakRound();
     666          if (done) break;
     667        }
     668
     669        // Saturate edges not satisfying the optimality condition
     670        Capacity delta;
     671        for (EdgeIt e(_graph); e != INVALID; ++e) {
     672          if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
     673            delta = _capacity[e] - (*_flow)[e];
     674            _excess[_graph.source(e)] -= delta;
     675            _excess[_graph.target(e)] += delta;
     676            (*_flow)[e] = _capacity[e];
     677          }
     678          if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
     679            _excess[_graph.target(e)] -= (*_flow)[e];
     680            _excess[_graph.source(e)] += (*_flow)[e];
     681            (*_flow)[e] = 0;
     682          }
     683        }
     684
     685        // Find active nodes (i.e. nodes with positive excess)
     686        for (NodeIt n(_graph); n != INVALID; ++n) {
     687          if (_excess[n] > 0) active_nodes.push_back(n);
     688        }
     689
     690        // Initialize the next edge maps
     691        for (NodeIt n(_graph); n != INVALID; ++n) {
     692          next_out[n] = OutEdgeIt(_graph, n);
     693          next_in[n] = InEdgeIt(_graph, n);
     694          next_dir[n] = true;
     695        }
     696
     697        // Perform push and relabel operations
     698        while (active_nodes.size() > 0) {
     699          // Select an active node (FIFO selection)
    517700          Node n = active_nodes[0], t;
    518701          bool relabel_enabled = true;
    519702
    520           // Performing push operations if there are admissible edges
    521           if (_excess[n] > 0) {
    522             for (OutEdgeIt e(_graph, n); e != INVALID; ++e) {
     703          // Perform push operations if there are admissible edges
     704          if (_excess[n] > 0 && next_dir[n]) {
     705            OutEdgeIt e = next_out[n];
     706            for ( ; e != INVALID; ++e) {
    523707              if (_capacity[e] - (*_flow)[e] > 0 && (*_red_cost)[e] < 0) {
    524                 delta = _capacity[e] - (*_flow)[e] <= _excess[n] ?
    525                         _capacity[e] - (*_flow)[e] : _excess[n];
     708                delta = std::min(_capacity[e] - (*_flow)[e], _excess[n]);
    526709                t = _graph.target(e);
    527710
     
    538721                if (ahead < 0) ahead = 0;
    539722
    540                 // Pushing flow along the edge
     723                // Push flow along the edge
    541724                if (ahead < delta) {
    542725                  (*_flow)[e] += ahead;
     
    558741              }
    559742            }
    560           }
    561 
    562           if (_excess[n] > 0) {
    563             for (InEdgeIt e(_graph, n); e != INVALID; ++e) {
     743            if (e != INVALID) {
     744              next_out[n] = e;
     745            } else {
     746              next_dir[n] = false;
     747            }
     748          }
     749
     750          if (_excess[n] > 0 && !next_dir[n]) {
     751            InEdgeIt e = next_in[n];
     752            for ( ; e != INVALID; ++e) {
    564753              if ((*_flow)[e] > 0 && -(*_red_cost)[e] < 0) {
    565                 delta = (*_flow)[e] <= _excess[n] ? (*_flow)[e] : _excess[n];
     754                delta = std::min((*_flow)[e], _excess[n]);
    566755                t = _graph.source(e);
    567756
     
    578767                if (ahead < 0) ahead = 0;
    579768
    580                 // Pushing flow along the edge
     769                // Push flow along the edge
    581770                if (ahead < delta) {
    582771                  (*_flow)[e] -= ahead;
     
    598787              }
    599788            }
    600           }
    601 
     789            next_in[n] = e;
     790          }
     791
     792          // Relabel the node if it is still active (or hyper)
    602793          if (relabel_enabled && (_excess[n] > 0 || hyper[n])) {
    603             // Performing relabel operation if the node is still active
    604             LCost min_red_cost = std::numeric_limits<LCost>::max();
     794            LCost min_red_cost = std::numeric_limits<LCost>::max() / 2;
    605795            for (OutEdgeIt oe(_graph, n); oe != INVALID; ++oe) {
    606796              if ( _capacity[oe] - (*_flow)[oe] > 0 &&
     
    614804            (*_potential)[n] -= min_red_cost + _epsilon;
    615805            hyper[n] = false;
    616           }
    617 
    618           // Removing active nodes with non-positive excess
     806
     807            // Reset the next edge maps
     808            next_out[n] = OutEdgeIt(_graph, n);
     809            next_in[n] = InEdgeIt(_graph, n);
     810            next_dir[n] = true;
     811          }
     812
     813          // Remove nodes that are not active nor hyper
    619814          while ( active_nodes.size() > 0 &&
    620815                  _excess[active_nodes[0]] <= 0 &&
     
    625820      }
    626821
    627       // Computing node potentials for the original costs
     822      // Compute node potentials for the original costs
    628823      ResidualCostMap<CostMap> res_cost(_orig_cost);
    629824      BellmanFord< ResGraph, ResidualCostMap<CostMap> >
     
    633828        (*_potential)[n] = bf.dist(n);
    634829
    635       // Handling non-zero lower bounds
     830      // Handle non-zero lower bounds
    636831      if (_lower) {
    637832        for (EdgeIt e(_graph); e != INVALID; ++e)
Note: See TracChangeset for help on using the changeset viewer.