RhodeCode

    /// Initializes the internal data structures and sets the initial

    /// flow to zero on each arc.

    void init() {

      createStructures();

      _phase = true;

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

        (*_excess)[n] = 0;

      }

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

        _flow->set(e, 0);

      }

      typename Digraph::template NodeMap<bool> reached(_graph, false);

      _level->initStart();

      _level->initAddItem(_target);

      std::vector<Node> queue;

      reached[_source] = true;

      queue.push_back(_target);

      reached[_target] = true;

      while (!queue.empty()) {

        _level->initNewLevel();

        std::vector<Node> nqueue;

        for (int i = 0; i < int(queue.size()); ++i) {

          Node n = queue[i];

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

            Node u = _graph.source(e);

            if (!reached[u] && _tolerance.positive((*_capacity)[e])) {

              reached[u] = true;

              _level->initAddItem(u);

              nqueue.push_back(u);

            }

          }

        }

        queue.swap(nqueue);

      }

      _level->initFinish();

      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {

        if (_tolerance.positive((*_capacity)[e])) {

          Node u = _graph.target(e);

          if ((*_level)[u] == _level->maxLevel()) continue;

          _flow->set(e, (*_capacity)[e]);

          (*_excess)[u] += (*_capacity)[e];

          if (u != _target && !_level->active(u)) {

            _level->activate(u);

          }

        }

      }

    }

    /// \brief Initializes the internal data structures using the

    /// given flow map.

    ///

    /// Initializes the internal data structures and sets the initial

    /// flow to the given \c flowMap. The \c flowMap should contain a

    /// flow or at least a preflow, i.e. at each node excluding the

    /// source node the incoming flow should greater or equal to the

    /// outgoing flow.

    /// \return \c false if the given \c flowMap is not a preflow.

    template <typename FlowMap>

    bool init(const FlowMap& flowMap) {

      createStructures();

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

        _flow->set(e, flowMap[e]);

      }

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

        Value excess = 0;

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

          excess += (*_flow)[e];

        }

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

          excess -= (*_flow)[e];

        }

        if (excess < 0 && n != _source) return false;

        (*_excess)[n] = excess;

      }

      typename Digraph::template NodeMap<bool> reached(_graph, false);

      _level->initStart();

      _level->initAddItem(_target);

      std::vector<Node> queue;

      reached[_source] = true;

      queue.push_back(_target);

      reached[_target] = true;

      while (!queue.empty()) {

        _level->initNewLevel();

        std::vector<Node> nqueue;

        for (int i = 0; i < int(queue.size()); ++i) {

          Node n = queue[i];

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

            Node u = _graph.source(e);

            if (!reached[u] &&

                _tolerance.positive((*_capacity)[e] - (*_flow)[e])) {

              reached[u] = true;

              _level->initAddItem(u);

              nqueue.push_back(u);

            }

          }

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

            Node v = _graph.target(e);

            if (!reached[v] && _tolerance.positive((*_flow)[e])) {

              reached[v] = true;

              _level->initAddItem(v);

              nqueue.push_back(v);

            }

          }

        }

        queue.swap(nqueue);

      }

      _level->initFinish();

      for (OutArcIt e(_graph, _source); e != INVALID; ++e) {

        Value rem = (*_capacity)[e] - (*_flow)[e];

        if (_tolerance.positive(rem)) {

          Node u = _graph.target(e);

          if ((*_level)[u] == _level->maxLevel()) continue;

          _flow->set(e, (*_capacity)[e]);

          (*_excess)[u] += rem;

        }

      }

      for (InArcIt e(_graph, _source); e != INVALID; ++e) {

        Value rem = (*_flow)[e];

        if (_tolerance.positive(rem)) {

          Node v = _graph.source(e);

          if ((*_level)[v] == _level->maxLevel()) continue;

          _flow->set(e, 0);

          (*_excess)[v] += rem;

        }

      }

      return true;

    }

    /// \brief Starts the first phase of the preflow algorithm.

    ///

    /// The preflow algorithm consists of two phases, this method runs

    /// the first phase. After the first phase the maximum flow value

    /// and a minimum value cut can already be computed, although a

    /// maximum flow is not yet obtained. So after calling this method

    /// \ref flowValue() returns the value of a maximum flow and \ref

    /// minCut() returns a minimum cut.

    /// \pre One of the \ref init() functions must be called before

    /// using this function.

    void startFirstPhase() {

      _phase = true;

      while (true) {

        int num = _node_num;

        Node n = INVALID;

        int level = -1;

        while (num > 0) {

          n = _level->highestActive();

          if (n == INVALID) goto first_phase_done;

          level = _level->highestActiveLevel();

          --num;

          Value excess = (*_excess)[n];

          int new_level = _level->maxLevel();

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

            Value rem = (*_capacity)[e] - (*_flow)[e];

            if (!_tolerance.positive(rem)) continue;

            Node v = _graph.target(e);

            if ((*_level)[v] < level) {

              if (!_level->active(v) && v != _target) {

                _level->activate(v);

              }

              if (!_tolerance.less(rem, excess)) {

                _flow->set(e, (*_flow)[e] + excess);

                (*_excess)[v] += excess;

                excess = 0;

                goto no_more_push_1;

              } else {

                excess -= rem;

                (*_excess)[v] += rem;

                _flow->set(e, (*_capacity)[e]);

              }

            } else if (new_level > (*_level)[v]) {

              new_level = (*_level)[v];

            }

          }

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

            Value rem = (*_flow)[e];

            if (!_tolerance.positive(rem)) continue;

            Node v = _graph.source(e);

            if ((*_level)[v] < level) {

              if (!_level->active(v) && v != _target) {

                _level->activate(v);

              }

              if (!_tolerance.less(rem, excess)) {

                _flow->set(e, (*_flow)[e] - excess);

                (*_excess)[v] += excess;

                excess = 0;

                goto no_more_push_1;

              } else {

                excess -= rem;

                (*_excess)[v] += rem;

                _flow->set(e, 0);

              }

            } else if (new_level > (*_level)[v]) {

              new_level = (*_level)[v];

            }

          }

        no_more_push_1:

          (*_excess)[n] = excess;

          if (excess != 0) {

            if (new_level + 1 < _level->maxLevel()) {

              _level->liftHighestActive(new_level + 1);

            } else {

              _level->liftHighestActiveToTop();

            }

            if (_level->emptyLevel(level)) {

              _level->liftToTop(level);

            }

          } else {

            _level->deactivate(n);

          }

        }

        num = _node_num * 20;

        while (num > 0) {

          while (level >= 0 && _level->activeFree(level)) {

            --level;

          }

          if (level == -1) {

            n = _level->highestActive();

            level = _level->highestActiveLevel();

            if (n == INVALID) goto first_phase_done;

          } else {

            n = _level->activeOn(level);

          }

          --num;

          Value excess = (*_excess)[n];

          int new_level = _level->maxLevel();

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

            Value rem = (*_capacity)[e] - (*_flow)[e];

            if (!_tolerance.positive(rem)) continue;

            Node v = _graph.target(e);

            if ((*_level)[v] < level) {

              if (!_level->active(v) && v != _target) {

                _level->activate(v);

              }

              if (!_tolerance.less(rem, excess)) {

                _flow->set(e, (*_flow)[e] + excess);

                (*_excess)[v] += excess;

                excess = 0;

                goto no_more_push_2;

              } else {

                excess -= rem;

                (*_excess)[v] += rem;

char test_lgf[] =

  "@nodes\n"

  "label\n"

  "0\n"

  "1\n"

  "2\n"

  "3\n"

  "4\n"

  "5\n"

  "6\n"

  "7\n"

  "8\n"

  "9\n"

  "@arcs\n"

  "    label capacity\n"

  "0 1 0     20\n"

  "0 2 1     0\n"

  "1 1 2     3\n"

  "1 2 3     8\n"

  "1 3 4     8\n"

  "2 5 5     5\n"

  "3 2 6     5\n"

  "3 5 7     5\n"

  "3 6 8     5\n"

  "4 3 9     3\n"

  "5 7 10    3\n"

  "5 6 11    10\n"

  "5 8 12    10\n"

  "6 8 13    8\n"

  "8 9 14    20\n"

  "8 1 15    5\n"

  "9 5 16    5\n"

  "@attributes\n"

  "source 1\n"

  "target 8\n";

void checkPreflowCompile()

{

  typedef int VType;

  typedef concepts::Digraph Digraph;

  typedef Digraph::Node Node;

  typedef Digraph::Arc Arc;

  typedef concepts::ReadMap<Arc,VType> CapMap;

  typedef concepts::ReadWriteMap<Arc,VType> FlowMap;

  typedef concepts::WriteMap<Node,bool> CutMap;

  typedef Elevator<Digraph, Digraph::Node> Elev;

  typedef LinkedElevator<Digraph, Digraph::Node> LinkedElev;

  Digraph g;

  Node n;

  Arc e;

  CapMap cap;

  FlowMap flow;

  CutMap cut;

  VType v;

  bool b;

  typedef Preflow<Digraph, CapMap>

            ::SetFlowMap<FlowMap>

            ::SetElevator<Elev>

            ::SetStandardElevator<LinkedElev>

            ::Create PreflowType;

  PreflowType preflow_test(g, cap, n, n);

  const PreflowType& const_preflow_test = preflow_test;

  const PreflowType::Elevator& elev = const_preflow_test.elevator();

  preflow_test.elevator(const_cast<PreflowType::Elevator&>(elev));

  PreflowType::Tolerance tol = const_preflow_test.tolerance();

  preflow_test.tolerance(tol);

  preflow_test

    .capacityMap(cap)

    .flowMap(flow)

    .source(n)

    .target(n);

  preflow_test.init();

  preflow_test.init(cap);

  preflow_test.startFirstPhase();

  preflow_test.startSecondPhase();

  preflow_test.run();

  preflow_test.runMinCut();

  v = const_preflow_test.flowValue();

  v = const_preflow_test.flow(e);

  const FlowMap& fm = const_preflow_test.flowMap();

  b = const_preflow_test.minCut(n);

  const_preflow_test.minCutMap(cut);

  ignore_unused_variable_warning(fm);

}

int cutValue (const SmartDigraph& g,

              const SmartDigraph::NodeMap<bool>& cut,

              const SmartDigraph::ArcMap<int>& cap) {

  int c=0;

  for(SmartDigraph::ArcIt e(g); e!=INVALID; ++e) {

    if (cut[g.source(e)] && !cut[g.target(e)]) c+=cap[e];

  }

  return c;

}

bool checkFlow(const SmartDigraph& g,

               const SmartDigraph::ArcMap<int>& flow,

               const SmartDigraph::ArcMap<int>& cap,

               SmartDigraph::Node s, SmartDigraph::Node t) {

  for (SmartDigraph::ArcIt e(g); e != INVALID; ++e) {

    if (flow[e] < 0 || flow[e] > cap[e]) return false;

  }

  for (SmartDigraph::NodeIt n(g); n != INVALID; ++n) {

    if (n == s || n == t) continue;

    int sum = 0;

    for (SmartDigraph::OutArcIt e(g, n); e != INVALID; ++e) {

      sum += flow[e];

    }

    for (SmartDigraph::InArcIt e(g, n); e != INVALID; ++e) {

      sum -= flow[e];

    }

    if (sum != 0) return false;

  }

  return true;

}

int main() {

  typedef SmartDigraph Digraph;

  typedef Digraph::Node Node;

  typedef Digraph::NodeIt NodeIt;

  typedef Digraph::ArcIt ArcIt;

  typedef Digraph::ArcMap<int> CapMap;

  typedef Digraph::ArcMap<int> FlowMap;

  typedef Digraph::NodeMap<bool> CutMap;

  typedef Preflow<Digraph, CapMap> PType;

  Digraph g;

  Node s, t;

  CapMap cap(g);

  std::istringstream input(test_lgf);

  DigraphReader<Digraph>(g,input).

    arcMap("capacity", cap).

    node("source",s).

    node("target",t).

    run();

  PType preflow_test(g, cap, s, t);

  preflow_test.run();

  check(checkFlow(g, preflow_test.flowMap(), cap, s, t),

        "The flow is not feasible.");

  CutMap min_cut(g);

  preflow_test.minCutMap(min_cut);

  int min_cut_value=cutValue(g,min_cut,cap);

  check(preflow_test.flowValue() == min_cut_value,

        "The max flow value is not equal to the three min cut values.");

  FlowMap flow(g);

  for(ArcIt e(g); e!=INVALID; ++e) flow[e] = preflow_test.flowMap()[e];

  int flow_value=preflow_test.flowValue();

  for(ArcIt e(g); e!=INVALID; ++e) cap[e]=2*cap[e];

  preflow_test.init(flow);

  preflow_test.startFirstPhase();

  CutMap min_cut1(g);

  preflow_test.minCutMap(min_cut1);

  min_cut_value=cutValue(g,min_cut1,cap);

  check(preflow_test.flowValue() == min_cut_value &&

        min_cut_value == 2*flow_value,

        "The max flow value or the min cut value is wrong.");

  preflow_test.startSecondPhase();

  check(checkFlow(g, preflow_test.flowMap(), cap, s, t),

        "The flow is not feasible.");

  CutMap min_cut2(g);

  preflow_test.minCutMap(min_cut2);

  min_cut_value=cutValue(g,min_cut2,cap);

  check(preflow_test.flowValue() == min_cut_value &&

        min_cut_value == 2*flow_value,

        "The max flow value or the three min cut values were not doubled");

  preflow_test.flowMap(flow);

  NodeIt tmp1(g,s);

  ++tmp1;

  if ( tmp1 != INVALID ) s=tmp1;

  NodeIt tmp2(g,t);

  ++tmp2;

  if ( tmp2 != INVALID ) t=tmp2;

  preflow_test.source(s);

  preflow_test.target(t);

  preflow_test.run();

  CutMap min_cut3(g);

  preflow_test.minCutMap(min_cut3);

  min_cut_value=cutValue(g,min_cut3,cap);

  check(preflow_test.flowValue() == min_cut_value,

        "The max flow value or the three min cut values are incorrect.");

  return 0;

}