RhodeCode

        _local_flow = true;

      }

      if (!_level) {

        _level = Traits::createElevator(_g, _node_num);

        _local_level = true;

      }

      if (!_excess) {

        _excess = new ExcessMap(_g);

      }

    }

    void destroyStructures() {

      if (_local_flow) {

        delete _flow;

      }

      if (_local_level) {

        delete _level;

      }

      if (_excess) {

        delete _excess;

      }

    }

  public:

    /// Sets the lower bound map.

    /// Sets the lower bound map.

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

    Circulation& lowerMap(const LowerMap& map) {

      _lo = ↦

      return *this;

    }

    /// Sets the upper bound (capacity) map.

    /// Sets the upper bound (capacity) map.

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

    Circulation& upperMap(const UpperMap& map) {

      _up = ↦

      return *this;

    }

    /// Sets the supply map.

    /// Sets the supply map.

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

    Circulation& supplyMap(const SupplyMap& map) {

      _supply = ↦

      return *this;

    }

    /// \brief Sets the flow map.

    ///

    /// Sets the flow map.

    /// If you don't use this function before calling \ref run() or

    /// \ref init(), an instance will be allocated automatically.

    /// The destructor deallocates this automatically allocated map,

    /// of course.

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

    Circulation& flowMap(FlowMap& map) {

      if (_local_flow) {

        delete _flow;

        _local_flow = false;

      }

      _flow = ↦

      return *this;

    }

    /// \brief Sets the elevator used by algorithm.

    ///

    /// Sets the elevator used by algorithm.

    /// If you don't use this function before calling \ref run() or

    /// \ref init(), an instance will be allocated automatically.

    /// The destructor deallocates this automatically allocated elevator,

    /// of course.

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

    Circulation& elevator(Elevator& elevator) {

      if (_local_level) {

        delete _level;

        _local_level = false;

      }

      _level = &elevator;

      return *this;

    }

    /// \brief Returns a const reference to the elevator.

    ///

    /// Returns a const reference to the elevator.

    ///

    /// \pre Either \ref run() or \ref init() must be called before

    /// using this function.

    const Elevator& elevator() const {

      return *_level;

    }

    ///

      _tol = tolerance;

      return *this;

    }

    /// \brief Returns a const reference to the tolerance.

    ///

    const Tolerance& tolerance() const {

    }

    /// \name Execution Control

    /// The simplest way to execute the algorithm is to call \ref run().\n

    /// If you need more control on the initial solution or the execution,

    /// first you have to call one of the \ref init() functions, then

    /// the \ref start() function.

    ///@{

    /// Initializes the internal data structures.

    /// Initializes the internal data structures and sets all flow values

    /// to the lower bound.

    void init()

    {

      LEMON_DEBUG(checkBoundMaps(),

        "Upper bounds must be greater or equal to the lower bounds");

      createStructures();

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

        (*_excess)[n] = (*_supply)[n];

      }

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

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

        (*_excess)[_g.target(e)] += (*_flow)[e];

        (*_excess)[_g.source(e)] -= (*_flow)[e];

      }

      // global relabeling tested, but in general case it provides

      // worse performance for random digraphs

      _level->initStart();

      for(NodeIt n(_g);n!=INVALID;++n)

        _level->initAddItem(n);

      _level->initFinish();

      for(NodeIt n(_g);n!=INVALID;++n)

        if(_tol.positive((*_excess)[n]))

          _level->activate(n);

    }

    /// Initializes the internal data structures using a greedy approach.

    /// Initializes the internal data structures using a greedy approach

    /// to construct the initial solution.

    void greedyInit()

    {

      LEMON_DEBUG(checkBoundMaps(),

        "Upper bounds must be greater or equal to the lower bounds");

      createStructures();

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

        (*_excess)[n] = (*_supply)[n];

      }

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

        if (!_tol.less(-(*_excess)[_g.target(e)], (*_up)[e])) {

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

          (*_excess)[_g.target(e)] += (*_up)[e];

          (*_excess)[_g.source(e)] -= (*_up)[e];

        } else if (_tol.less(-(*_excess)[_g.target(e)], (*_lo)[e])) {

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

          (*_excess)[_g.target(e)] += (*_lo)[e];

          (*_excess)[_g.source(e)] -= (*_lo)[e];

        } else {

          Value fc = -(*_excess)[_g.target(e)];

          _flow->set(e, fc);

          (*_excess)[_g.target(e)] = 0;

          (*_excess)[_g.source(e)] -= fc;

        }

      }

      _level->initStart();

      for(NodeIt n(_g);n!=INVALID;++n)

        _level->initAddItem(n);

      _level->initFinish();

      for(NodeIt n(_g);n!=INVALID;++n)

        if(_tol.positive((*_excess)[n]))

          _level->activate(n);

    }

    ///Executes the algorithm

    ///This function executes the algorithm.

    ///

    ///\return \c true if a feasible circulation is found.

    ///

    ///\sa barrier()

    ///\sa barrierMap()

    bool start()

    {

      Node act;

      Node bact=INVALID;

 *

 * Copyright (C) 2003-2009

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#ifndef LEMON_PREFLOW_H

#define LEMON_PREFLOW_H

#include <lemon/tolerance.h>

#include <lemon/elevator.h>

/// \file

/// \ingroup max_flow

/// \brief Implementation of the preflow algorithm.

namespace lemon {

  /// \brief Default traits class of Preflow class.

  ///

  /// Default traits class of Preflow class.

  /// \tparam GR Digraph type.

  /// \tparam CAP Capacity map type.

  template <typename GR, typename CAP>

  struct PreflowDefaultTraits {

    /// \brief The type of the digraph the algorithm runs on.

    typedef GR Digraph;

    /// \brief The type of the map that stores the arc capacities.

    ///

    /// The type of the map that stores the arc capacities.

    /// It must meet the \ref concepts::ReadMap "ReadMap" concept.

    typedef CAP CapacityMap;

    /// \brief The type of the flow values.

    typedef typename CapacityMap::Value Value;

    /// \brief The type of the map that stores the flow values.

    ///

    /// The type of the map that stores the flow values.

    /// It must meet the \ref concepts::ReadWriteMap "ReadWriteMap" concept.

    typedef typename Digraph::template ArcMap<Value> FlowMap;

    /// \brief Instantiates a FlowMap.

    ///

    /// This function instantiates a \ref FlowMap.

    /// \param digraph The digraph for which we would like to define

    /// the flow map.

    static FlowMap* createFlowMap(const Digraph& digraph) {

      return new FlowMap(digraph);

    }

    /// \brief The elevator type used by Preflow algorithm.

    ///

    /// The elevator type used by Preflow algorithm.

    ///

    /// \sa Elevator

    /// \sa LinkedElevator

    typedef LinkedElevator<Digraph, typename Digraph::Node> Elevator;

    /// \brief Instantiates an Elevator.

    ///

    /// This function instantiates an \ref Elevator.

    /// \param digraph The digraph for which we would like to define

    /// the elevator.

    /// \param max_level The maximum level of the elevator.

    static Elevator* createElevator(const Digraph& digraph, int max_level) {

      return new Elevator(digraph, max_level);

    }

    /// \brief The tolerance used by the algorithm

    ///

    /// The tolerance used by the algorithm to handle inexact computation.

    typedef lemon::Tolerance<Value> Tolerance;

  };

  /// \ingroup max_flow

  ///

  /// \brief %Preflow algorithm class.

  ///

  /// This class provides an implementation of Goldberg-Tarjan's \e preflow

  /// \e push-relabel algorithm producing a \ref max_flow

  /// "flow of maximum value" in a digraph.

  /// The preflow algorithms are the fastest known maximum

  /// \e "highest label" and the \e "bound decrease" heuristics.

  /// The worst case time complexity of the algorithm is \f$O(n^2\sqrt{e})\f$.

  ///

  /// The algorithm consists of two phases. After the first phase

  /// the maximum flow value and the minimum cut is obtained. The

  /// second phase constructs a feasible maximum flow on each arc.

  ///

  /// \tparam GR The type of the digraph the algorithm runs on.

  /// \tparam CAP The type of the capacity map. The default map

  /// type is \ref concepts::Digraph::ArcMap "GR::ArcMap<int>".

#ifdef DOXYGEN

  template <typename GR, typename CAP, typename TR>

#else

  template <typename GR,

            typename CAP = typename GR::template ArcMap<int>,

            typename TR = PreflowDefaultTraits<GR, CAP> >

#endif

  class Preflow {

  public:

    ///The \ref PreflowDefaultTraits "traits class" of the algorithm.

    typedef TR Traits;

    ///The type of the digraph the algorithm runs on.

    typedef typename Traits::Digraph Digraph;

    ///The type of the capacity map.

    typedef typename Traits::CapacityMap CapacityMap;

    ///The type of the flow values.

    typedef typename Traits::Value Value;

    ///The type of the flow map.

    typedef typename Traits::FlowMap FlowMap;

    ///The type of the elevator.

    typedef typename Traits::Elevator Elevator;

    ///The type of the tolerance.

    typedef typename Traits::Tolerance Tolerance;

  private:

    TEMPLATE_DIGRAPH_TYPEDEFS(Digraph);

    const Digraph& _graph;

    const CapacityMap* _capacity;

    int _node_num;

    Node _source, _target;

    FlowMap* _flow;

    bool _local_flow;

    Elevator* _level;

    bool _local_level;

    typedef typename Digraph::template NodeMap<Value> ExcessMap;

    ExcessMap* _excess;

    Tolerance _tolerance;

    bool _phase;

    void createStructures() {

      _node_num = countNodes(_graph);

      if (!_flow) {

        _flow = Traits::createFlowMap(_graph);

        _local_flow = true;

      }

      if (!_level) {

        _level = Traits::createElevator(_graph, _node_num);

        _local_level = true;

      }

      if (!_excess) {

        _excess = new ExcessMap(_graph);

      }

    }

    void destroyStructures() {

      if (_local_flow) {

        delete _flow;

      }

      if (_local_level) {

        delete _level;

      }

      if (_excess) {

        delete _excess;

      }

    }

  public:

    typedef Preflow Create;

    ///\name Named Template Parameters

    ///@{

  public:

    /// \brief The constructor of the class.

    ///

    /// The constructor of the class.

    /// \param digraph The digraph the algorithm runs on.

    /// \param capacity The capacity of the arcs.

    /// \param source The source node.

    /// \param target The target node.

    Preflow(const Digraph& digraph, const CapacityMap& capacity,

            Node source, Node target)

      : _graph(digraph), _capacity(&capacity),

        _node_num(0), _source(source), _target(target),

        _flow(0), _local_flow(false),

        _level(0), _local_level(false),

        _excess(0), _tolerance(), _phase() {}

    /// \brief Destructor.

    ///

    /// Destructor.

    ~Preflow() {

      destroyStructures();

    }

    /// \brief Sets the capacity map.

    ///

    /// Sets the capacity map.

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

    Preflow& capacityMap(const CapacityMap& map) {

      _capacity = ↦

      return *this;

    }

    /// \brief Sets the flow map.

    ///

    /// Sets the flow map.

    /// If you don't use this function before calling \ref run() or

    /// \ref init(), an instance will be allocated automatically.

    /// The destructor deallocates this automatically allocated map,

    /// of course.

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

    Preflow& flowMap(FlowMap& map) {

      if (_local_flow) {

        delete _flow;

        _local_flow = false;

      }

      _flow = ↦

      return *this;

    }

    /// \brief Sets the source node.

    ///

    /// Sets the source node.

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

    Preflow& source(const Node& node) {

      _source = node;

      return *this;

    }

    /// \brief Sets the target node.

    ///

    /// Sets the target node.

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

    Preflow& target(const Node& node) {

      _target = node;

      return *this;

    }

    /// \brief Sets the elevator used by algorithm.

    ///

    /// Sets the elevator used by algorithm.

    /// If you don't use this function before calling \ref run() or

    /// \ref init(), an instance will be allocated automatically.

    /// The destructor deallocates this automatically allocated elevator,

    /// of course.

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

    Preflow& elevator(Elevator& elevator) {

      if (_local_level) {

        delete _level;

        _local_level = false;

      }

      _level = &elevator;

      return *this;

    }

    /// \brief Returns a const reference to the elevator.

    ///

    /// Returns a const reference to the elevator.

    ///

    /// \pre Either \ref run() or \ref init() must be called before

    /// using this function.

    const Elevator& elevator() const {

      return *_level;

    }

    ///

      _tolerance = tolerance;

      return *this;

    }

    /// \brief Returns a const reference to the tolerance.

    ///

    const Tolerance& tolerance() const {

    }

    /// \name Execution Control

    /// The simplest way to execute the preflow algorithm is to use

    /// \ref run() or \ref runMinCut().\n

    /// If you need more control on the initial solution or the execution,

    /// first you have to call one of the \ref init() functions, then

    /// \ref startFirstPhase() and if you need it \ref startSecondPhase().

    ///@{

    /// \brief Initializes the internal data structures.

    ///

    /// 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;

      }

/* -*- mode: C++; indent-tabs-mode: nil; -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library.

 *

 * Copyright (C) 2003-2009

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#include <iostream>

#include "test_tools.h"

#include <lemon/list_graph.h>

#include <lemon/circulation.h>

#include <lemon/lgf_reader.h>

#include <lemon/concepts/digraph.h>

#include <lemon/concepts/maps.h>

using namespace lemon;

char test_lgf[] =

  "@nodes\n"

  "label\n"

  "0\n"

  "1\n"

  "2\n"

  "3\n"

  "4\n"

  "5\n"

  "@arcs\n"

  "     lcap  ucap\n"

  "0 1  2  10\n"

  "0 2  2  6\n"

  "1 3  4  7\n"

  "1 4  0  5\n"

  "2 4  1  3\n"

  "3 5  3  8\n"

  "4 5  3  7\n"

  "@attributes\n"

  "source 0\n"

  "sink   5\n";

void checkCirculationCompile()

{

  typedef int VType;

  typedef concepts::Digraph Digraph;

  typedef Digraph::Node Node;

  typedef Digraph::Arc Arc;

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

  typedef concepts::ReadMap<Node,VType> SupplyMap;

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

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

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

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

  Digraph g;

  Node n;

  Arc a;

  CapMap lcap, ucap;

  SupplyMap supply;

  FlowMap flow;

  BarrierMap bar;

  VType v;

  bool b;

  typedef Circulation<Digraph, CapMap, CapMap, SupplyMap>

            ::SetFlowMap<FlowMap>

            ::SetElevator<Elev>

            ::SetStandardElevator<LinkedElev>

            ::Create CirculationType;

  CirculationType circ_test(g, lcap, ucap, supply);

  const CirculationType& const_circ_test = circ_test;

  circ_test

    .lowerMap(lcap)

    .upperMap(ucap)

    .supplyMap(supply)

    .flowMap(flow);

  circ_test.init();

  circ_test.greedyInit();

  circ_test.start();

  circ_test.run();

  v = const_circ_test.flow(a);

  const FlowMap& fm = const_circ_test.flowMap();

  b = const_circ_test.barrier(n);

  const_circ_test.barrierMap(bar);

  ignore_unused_variable_warning(fm);

}

template <class G, class LM, class UM, class DM>

void checkCirculation(const G& g, const LM& lm, const UM& um,

                      const DM& dm, bool find)

{

  Circulation<G, LM, UM, DM> circ(g, lm, um, dm);

  bool ret = circ.run();

  if (find) {

    check(ret, "A feasible solution should have been found.");

    check(circ.checkFlow(), "The found flow is corrupt.");

    check(!circ.checkBarrier(), "A barrier should not have been found.");

  } else {

    check(!ret, "A feasible solution should not have been found.");

    check(circ.checkBarrier(), "The found barrier is corrupt.");

  }

}

int main (int, char*[])

{

  typedef ListDigraph Digraph;

  DIGRAPH_TYPEDEFS(Digraph);

  Digraph g;

  IntArcMap lo(g), up(g);

  IntNodeMap delta(g, 0);

  Node s, t;

  std::istringstream input(test_lgf);

  DigraphReader<Digraph>(g,input).

    arcMap("lcap", lo).

    arcMap("ucap", up).

    node("source",s).

    node("sink",t).

    run();

  delta[s] = 7; delta[t] = -7;

  checkCirculation(g, lo, up, delta, true);

  delta[s] = 13; delta[t] = -13;

  checkCirculation(g, lo, up, delta, true);

  delta[s] = 6; delta[t] = -6;

  checkCirculation(g, lo, up, delta, false);

  delta[s] = 14; delta[t] = -14;

  checkCirculation(g, lo, up, delta, false);

  delta[s] = 7; delta[t] = -13;

  checkCirculation(g, lo, up, delta, true);

  delta[s] = 5; delta[t] = -15;

  checkCirculation(g, lo, up, delta, true);

  delta[s] = 10; delta[t] = -11;

  checkCirculation(g, lo, up, delta, true);

  delta[s] = 11; delta[t] = -10;

  checkCirculation(g, lo, up, delta, false);

  return 0;

}

/* -*- mode: C++; indent-tabs-mode: nil; -*-

 *

 * This file is a part of LEMON, a generic C++ optimization library.

 *

 * Copyright (C) 2003-2009

 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport

 * (Egervary Research Group on Combinatorial Optimization, EGRES).

 *

 * Permission to use, modify and distribute this software is granted

 * provided that this copyright notice appears in all copies. For

 * precise terms see the accompanying LICENSE file.

 *

 * This software is provided "AS IS" with no warranty of any kind,

 * express or implied, and with no claim as to its suitability for any

 * purpose.

 *

 */

#include <iostream>

#include "test_tools.h"

#include <lemon/smart_graph.h>

#include <lemon/preflow.h>

#include <lemon/concepts/digraph.h>

#include <lemon/concepts/maps.h>

#include <lemon/lgf_reader.h>

#include <lemon/elevator.h>

using namespace lemon;

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;

  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];