RhodeCode

    /// \ref resetParams() or \ref reset() is used.

    /// If the underlying digraph was also modified after the construction

    /// of the class or the last \ref reset() call, then the \ref reset()

    /// function must be used, otherwise \ref resetParams() is sufficient.

    ///

    /// For example,

    /// \code

    ///   CostScaling<ListDigraph> cs(graph);

    ///

    ///   // First run

    ///   cs.lowerMap(lower).upperMap(upper).costMap(cost)

    ///     .supplyMap(sup).run();

    ///

    ///   // Run again with modified cost map (resetParams() is not called,

    ///   // so only the cost map have to be set again)

    ///   cost[e] += 100;

    ///   cs.costMap(cost).run();

    ///

    ///   // Run again from scratch using resetParams()

    ///   // (the lower bounds will be set to zero on all arcs)

    ///   cs.resetParams();

    ///   cs.upperMap(capacity).costMap(cost)

    ///     .supplyMap(sup).run();

    /// \endcode

    ///

    /// \return <tt>(*this)</tt>

    ///

    /// \see reset(), run()

    CostScaling& resetParams() {

      for (int i = 0; i != _res_node_num; ++i) {

        _supply[i] = 0;

      }

      int limit = _first_out[_root];

      for (int j = 0; j != limit; ++j) {

        _lower[j] = 0;

        _upper[j] = INF;

        _scost[j] = _forward[j] ? 1 : -1;

      }

      for (int j = limit; j != _res_arc_num; ++j) {

        _lower[j] = 0;

        _upper[j] = INF;

        _scost[j] = 0;

        _scost[_reverse[j]] = 0;

      }

      _have_lower = false;

      return *this;

    }

    ///

    ///

    /// \return <tt>(*this)</tt>

    CostScaling& reset() {

      // Resize vectors

      _node_num = countNodes(_graph);

      _arc_num = countArcs(_graph);

      _res_node_num = _node_num + 1;

      _res_arc_num = 2 * (_arc_num + _node_num);

      _root = _node_num;

      _first_out.resize(_res_node_num + 1);

      _forward.resize(_res_arc_num);

      _source.resize(_res_arc_num);

      _target.resize(_res_arc_num);

      _reverse.resize(_res_arc_num);

      _lower.resize(_res_arc_num);

      _upper.resize(_res_arc_num);

      _scost.resize(_res_arc_num);

      _supply.resize(_res_node_num);

      _res_cap.resize(_res_arc_num);

      _cost.resize(_res_arc_num);

      _pi.resize(_res_node_num);

      _excess.resize(_res_node_num);

      _next_out.resize(_res_node_num);

      _arc_vec.reserve(_res_arc_num);

      _cost_vec.reserve(_res_arc_num);

      // Copy the graph

      int i = 0, j = 0, k = 2 * _arc_num + _node_num;

      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {

        _node_id[n] = i;

      }

      i = 0;

      for (NodeIt n(_graph); n != INVALID; ++n, ++i) {

        _first_out[i] = j;

        for (OutArcIt a(_graph, n); a != INVALID; ++a, ++j) {

          _arc_idf[a] = j;

          _forward[j] = true;

          _source[j] = i;

          _target[j] = _node_id[_graph.runningNode(a)];

        }

        for (InArcIt a(_graph, n); a != INVALID; ++a, ++j) {

          _arc_idb[a] = j;

          _forward[j] = false;

          _source[j] = i;

          _target[j] = _node_id[_graph.runningNode(a)];

        }

      }

      _sup_node_num = 0;

      for (NodeIt n(_graph); n != INVALID; ++n) {

        if (sup[n] > 0) ++_sup_node_num;

      }

      // Find a feasible flow using Circulation

      Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>

        circ(_graph, low, cap, sup);

      if (!circ.flowMap(flow).run()) return INFEASIBLE;

      // Set residual capacities and handle GEQ supply type

      if (_sum_supply < 0) {

        for (ArcIt a(_graph); a != INVALID; ++a) {

          Value fa = flow[a];

          _res_cap[_arc_idf[a]] = cap[a] - fa;

          _res_cap[_arc_idb[a]] = fa;

          sup[_graph.source(a)] -= fa;

          sup[_graph.target(a)] += fa;

        }

        for (NodeIt n(_graph); n != INVALID; ++n) {

          _excess[_node_id[n]] = sup[n];

        }

        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {

          int u = _target[a];

          int ra = _reverse[a];

          _res_cap[a] = -_sum_supply + 1;

          _res_cap[ra] = -_excess[u];

          _cost[a] = 0;

          _cost[ra] = 0;

          _excess[u] = 0;

        }

      } else {

        for (ArcIt a(_graph); a != INVALID; ++a) {

          Value fa = flow[a];

          _res_cap[_arc_idf[a]] = cap[a] - fa;

          _res_cap[_arc_idb[a]] = fa;

        }

        for (int a = _first_out[_root]; a != _res_arc_num; ++a) {

          int ra = _reverse[a];

          _res_cap[a] = 0;

          _res_cap[ra] = 0;

          _cost[a] = 0;

          _cost[ra] = 0;

        }

      }

      // Initialize data structures for buckets

      _max_rank = _alpha * _res_node_num;

      _buckets.resize(_max_rank);

      _bucket_next.resize(_res_node_num + 1);

      _bucket_prev.resize(_res_node_num + 1);

      _rank.resize(_res_node_num + 1);

      switch (method) {

        case PUSH:

          startPush();

          break;

        case AUGMENT:

          startAugment(_res_node_num - 1);

          break;

        case PARTIAL_AUGMENT:

          break;

      }

      // Compute node potentials for the original costs

      _arc_vec.clear();

      _cost_vec.clear();

      for (int j = 0; j != _res_arc_num; ++j) {

        if (_res_cap[j] > 0) {

          _arc_vec.push_back(IntPair(_source[j], _target[j]));

          _cost_vec.push_back(_scost[j]);

        }

      }

      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());

      typename BellmanFord<StaticDigraph, LargeCostArcMap>

        ::template SetDistMap<LargeCostNodeMap>::Create bf(_sgr, _cost_map);

      bf.distMap(_pi_map);

      bf.init(0);

      bf.start();

      // Handle non-zero lower bounds

      if (_have_lower) {

        int limit = _first_out[_root];

        for (int j = 0; j != limit; ++j) {

          if (!_forward[j]) _res_cap[j] += _lower[j];

        }

      }

    }

    // Initialize a cost scaling phase

    void initPhase() {

      // Saturate arcs not satisfying the optimality condition

      for (int u = 0; u != _res_node_num; ++u) {

        int last_out = _first_out[u+1];

        LargeCost pi_u = _pi[u];

        for (int a = _first_out[u]; a != last_out; ++a) {

          int v = _target[a];

            _excess[u] -= delta;

            _excess[v] += delta;

            _res_cap[a] = 0;

            _res_cap[_reverse[a]] += delta;

          }

        }

      }

      // Find active nodes (i.e. nodes with positive excess)

      for (int u = 0; u != _res_node_num; ++u) {

        if (_excess[u] > 0) _active_nodes.push_back(u);

      }

      // Initialize the next arcs

      for (int u = 0; u != _res_node_num; ++u) {

        _next_out[u] = _first_out[u];

      }

    }

    // Early termination heuristic

    bool earlyTermination() {

      const double EARLY_TERM_FACTOR = 3.0;

      // Build a static residual graph

      _arc_vec.clear();

      _cost_vec.clear();

      for (int j = 0; j != _res_arc_num; ++j) {

        if (_res_cap[j] > 0) {

          _arc_vec.push_back(IntPair(_source[j], _target[j]));

          _cost_vec.push_back(_cost[j] + 1);

        }

      }

      _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());

      // Run Bellman-Ford algorithm to check if the current flow is optimal

      BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);

      bf.init(0);

      bool done = false;

      int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));

      for (int i = 0; i < K && !done; ++i) {

        done = bf.processNextWeakRound();

      }

      return done;

    }

    // Global potential update heuristic

    void globalUpdate() {

      // Initialize buckets

      for (int r = 0; r != _max_rank; ++r) {

        _buckets[r] = bucket_end;

      }

      Value total_excess = 0;

      for (int i = 0; i != _res_node_num; ++i) {

        if (_excess[i] < 0) {

          _rank[i] = 0;

        } else {

          total_excess += _excess[i];

          _rank[i] = _max_rank;

        }

      }

      if (total_excess == 0) return;

      // Search the buckets

      int r = 0;

      for ( ; r != _max_rank; ++r) {

        while (_buckets[r] != bucket_end) {

          // Remove the first node from the current bucket

          int u = _buckets[r];

          _buckets[r] = _bucket_next[u];

          // Search the incomming arcs of u

          LargeCost pi_u = _pi[u];

          int last_out = _first_out[u+1];

          for (int a = _first_out[u]; a != last_out; ++a) {

            int ra = _reverse[a];

            if (_res_cap[ra] > 0) {

              int v = _source[ra];

              int old_rank_v = _rank[v];

              if (r < old_rank_v) {

                // Compute the new rank of v

                LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;

                int new_rank_v = old_rank_v;

                // Change the rank of v

                if (new_rank_v < old_rank_v) {

                  _rank[v] = new_rank_v;

                  _next_out[v] = _first_out[v];

                  // Remove v from its old bucket

                  if (old_rank_v < _max_rank) {

                    if (_buckets[old_rank_v] == v) {

                      _buckets[old_rank_v] = _bucket_next[v];

                    } else {

                    }

                  }

                  _buckets[new_rank_v] = v;

                }

              }

            }

          }

          // Finish search if there are no more active nodes

          if (_excess[u] > 0) {

            total_excess -= _excess[u];

            if (total_excess <= 0) break;

          }

        }

        if (total_excess <= 0) break;

      }

      // Relabel nodes

      for (int u = 0; u != _res_node_num; ++u) {

        int k = std::min(_rank[u], r);

        if (k > 0) {

          _pi[u] -= _epsilon * k;

          _next_out[u] = _first_out[u];

        }

      }

    }

    /// Execute the algorithm performing augment and relabel operations

    void startAugment(int max_length) {

      // Paramters for heuristics

      const int EARLY_TERM_EPSILON_LIMIT = 1000;

      const double GLOBAL_UPDATE_FACTOR = 3.0;

      const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *

        (_res_node_num + _sup_node_num * _sup_node_num));

      int next_update_limit = global_update_freq;

      int relabel_cnt = 0;

      // Perform cost scaling phases

      std::vector<int> path;

      for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?

                                        1 : _epsilon / _alpha )

      {

        // Early termination heuristic

        if (_epsilon <= EARLY_TERM_EPSILON_LIMIT) {

          if (earlyTermination()) break;

        }

        // Initialize current phase