RhodeCode

digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and

A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the

\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]

\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)

    \quad \forall u\in V\setminus\{s,t\} \f]

\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]

This group contains the algorithms for finding minimum cost flows and

in a network with capacity constraints (lower and upper bounds)

and arc costs.

upper bounds for the flow values on the arcs, for which

\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.

on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the

signed supply values of the nodes.

If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$

supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with

\f$-sup(u)\f$ demand.

A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution

of the following optimization problem.

\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]

\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq

    sup(u) \quad \forall u\in V \f]

\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]

The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be

zero or negative in order to have a feasible solution (since the sum

of the expressions on the left-hand side of the inequalities is zero).

It means that the total demand must be greater or equal to the total

supply and all the supplies have to be carried out from the supply nodes,

but there could be demands that are not satisfied.

If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand

constraints have to be satisfied with equality, i.e. all demands

have to be satisfied and all supplies have to be used.

If you need the opposite inequalities in the supply/demand constraints

(i.e. the total demand is less than the total supply and all the demands

have to be satisfied while there could be supplies that are not used),

then you could easily transform the problem to the above form by reversing

the direction of the arcs and taking the negative of the supply values

(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).

However \ref NetworkSimplex algorithm also supports this form directly

for the sake of convenience.

A feasible solution for this problem can be found using \ref Circulation.

Note that the above formulation is actually more general than the usual

definition of the minimum cost flow problem, in which strict equalities

are required in the supply/demand contraints, i.e.

\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =

    sup(u) \quad \forall u\in V. \f]

However if the sum of the supply values is zero, then these two problems

are equivalent. So if you need the equality form, you have to ensure this

additional contraint for the algorithms.

The dual solution of the minimum cost flow problem is represented by node 

potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.

An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem

is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$

node potentials the following \e complementary \e slackness optimality

conditions hold.

 - For all \f$uv\in A\f$ arcs:

   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;

   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;

   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.

 - For all \f$u\in V\f$:

   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,

     then \f$\pi(u)=0\f$.

Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc

\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.

\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]

All algorithms provide dual solution (node potentials) as well

if an optimal flow is found.

LEMON contains several algorithms for solving minimum cost flow problems.

 - \ref NetworkSimplex Primal Network Simplex algorithm with various

   pivot strategies.

 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on

   cost scaling.

 - \ref CapacityScaling Successive Shortest %Path algorithm with optional

 - \ref CancelAndTighten The Cancel and Tighten algorithm.

 - \ref CycleCanceling Cycle-Canceling algorithms.

Most of these implementations support the general inequality form of the

minimum cost flow problem, but CancelAndTighten and CycleCanceling

only support the equality form due to the primal method they use.

In general NetworkSimplex is the most efficient implementation,

but in special cases other algorithms could be faster.

For example, if the total supply and/or capacities are rather small,

CapacityScaling is usually the fastest algorithm (without effective scaling).

#include <lemon/maps.h>

#include <lemon/circulation.h>

#include <lemon/adaptors.h>

  /// Moreover it supports both direction of the supply/demand inequality

  /// constraints. For more information see \ref ProblemType.

  /// implementations, from which the most efficient one is used

  /// by default. For more information see \ref PivotRule.

#ifdef DOXYGEN

    /// The type of the flow map

    typedef GR::ArcMap<Flow> FlowMap;

    /// The type of the potential map

    typedef GR::NodeMap<Cost> PotentialMap;

#else

#endif

    /// \brief Enum type for selecting the problem type.

    ///

    /// Enum type for selecting the problem type, i.e. the direction of

    /// the inequalities in the supply/demand constraints of the

    /// \ref min_cost_flow "minimum cost flow problem".

    ///

    /// The default problem type is \c GEQ, since this form is supported

    /// by other minimum cost flow algorithms and the \ref Circulation

    /// algorithm as well.

    /// The \c LEQ problem type can be selected using the \ref problemType()

    /// function.

    ///

    /// Note that the equality form is a special case of both problem type.

    enum ProblemType {

      /// This option means that there are "<em>greater or equal</em>"

      /// constraints in the defintion, i.e. the exact formulation of the

      /// problem is the following.

      /**

          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]

          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq

              sup(u) \quad \forall u\in V \f]

          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]

      */

      /// It means that the total demand must be greater or equal to the 

      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or

      /// negative) and all the supplies have to be carried out from 

      /// the supply nodes, but there could be demands that are not 

      /// satisfied.

      GEQ,

      /// It is just an alias for the \c GEQ option.

      CARRY_SUPPLIES = GEQ,

      /// This option means that there are "<em>less or equal</em>"

      /// constraints in the defintion, i.e. the exact formulation of the

      /// problem is the following.

      /**

          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]

          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq

              sup(u) \quad \forall u\in V \f]

          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]

      */

      /// It means that the total demand must be less or equal to the 

      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or

      /// positive) and all the demands have to be satisfied, but there

      /// could be supplies that are not carried out from the supply

      /// nodes.

      LEQ,

      /// It is just an alias for the \c LEQ option.

      SATISFY_DEMANDS = LEQ

    };

    ProblemType _ptype;

    /// The constructor of the class.

      _psupply(NULL), _pstsup(false), _ptype(GEQ),

    /// \name Parameters

    /// The parameters of the algorithm can be specified using these

    /// functions.

    /// @{

    /// \brief Set the problem type.

    ///

    /// This function sets the problem type for the algorithm.

    /// If it is not used before calling \ref run(), the \ref GEQ problem

    /// type will be used.

    ///

    /// For more information see \ref ProblemType.

    ///

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

    NetworkSimplex& problemType(ProblemType problem_type) {

      _ptype = problem_type;

      return *this;

    }

    /// @}

    /// The paramters can be specified using functions \ref lowerMap(),

    /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 

    /// \ref problemType(), \ref flowMap() and \ref potentialMap().

    /// For example,

    /// before using functions \ref lowerMap(), \ref upperMap(),

    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),

    /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 

    /// \ref flowMap() and \ref potentialMap().

      _ptype = GEQ;

      if (_local_flow) delete _flow_map;

      if (_local_potential) delete _potential_map;

      _flow_map = NULL;

      _potential_map = NULL;

      _local_flow = false;

      _local_potential = false;

      Flow sum_supply = 0;

          sum_supply += _supply[i];

        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||

                       (_ptype == LEQ && sum_supply >= 0);

        _supply[_node_id[_ptarget]] = -_pstflow;

      // Infinite capacity value

      Flow inf_cap =

        std::numeric_limits<Flow>::has_infinity ?

        std::numeric_limits<Flow>::infinity() :

        std::numeric_limits<Flow>::max();

      // Initialize artifical cost

      Cost art_cost;

      if (std::numeric_limits<Cost>::is_exact) {

        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;

      } else {

        art_cost = std::numeric_limits<Cost>::min();

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

          if (_cost[i] > art_cost) art_cost = _cost[i];

        }

        art_cost = (art_cost + 1) * _node_num;

      }

      // Run Circulation to check if a feasible solution exists

      typedef ConstMap<Arc, Flow> ConstArcMap;

      FlowNodeMap *csup = NULL;

      bool local_csup = false;

      if (_psupply) {

        csup = _psupply;

      } else {

        csup = new FlowNodeMap(_graph, 0);

        (*csup)[_psource] =  _pstflow;

        (*csup)[_ptarget] = -_pstflow;

        local_csup = true;

      }

      bool circ_result = false;

      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {

        // GEQ problem type

        if (_plower) {

          if (_pupper) {

            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>

              circ(_graph, *_plower, *_pupper, *csup);

            circ_result = circ.run();

          } else {

            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>

              circ(_graph, *_plower, ConstArcMap(inf_cap), *csup);

            circ_result = circ.run();

          }

        } else {

          if (_pupper) {

            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>

              circ(_graph, ConstArcMap(0), *_pupper, *csup);

            circ_result = circ.run();

          } else {

            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>

              circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *csup);

            circ_result = circ.run();

          }

        }

      } else {

        // LEQ problem type

        typedef ReverseDigraph<const GR> RevGraph;

        typedef NegMap<FlowNodeMap> NegNodeMap;

        RevGraph rgraph(_graph);

        NegNodeMap neg_csup(*csup);

        if (_plower) {

          if (_pupper) {

            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>

              circ(rgraph, *_plower, *_pupper, neg_csup);

            circ_result = circ.run();

          } else {

            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>

              circ(rgraph, *_plower, ConstArcMap(inf_cap), neg_csup);

            circ_result = circ.run();

          }

        } else {

          if (_pupper) {

            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>

              circ(rgraph, ConstArcMap(0), *_pupper, neg_csup);

            circ_result = circ.run();

          } else {

            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>

              circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup);

            circ_result = circ.run();

          }

        }

      }

      if (local_csup) delete csup;

      if (!circ_result) return false;

      _supply[_root] = -sum_supply;

      if (sum_supply < 0) {

        _pi[_root] = -art_cost;

      } else {

        _pi[_root] = art_cost;

      }

        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {

          _pi[u] = -art_cost + _pi[_root];

          _pi[u] = art_cost + _pi[_root];

  "label  sup1 sup2 sup3 sup4 sup5\n"

  "    1    20   27    0   20   30\n"

  "    2    -4    0    0   -8   -3\n"

  "    3     0    0    0    0    0\n"

  "    4     0    0    0    0    0\n"

  "    5     9    0    0    6   11\n"

  "    6    -6    0    0   -5   -6\n"

  "    7     0    0    0    0    0\n"

  "    8     0    0    0    0    3\n"

  "    9     3    0    0    0    0\n"

  "   10    -2    0    0   -7   -2\n"

  "   11     0    0    0  -10    0\n"

  "   12   -20  -27    0  -30  -20\n"

enum ProblemType {

  EQ,

  GEQ,

  LEQ

};

             .flowMap(flow)

             .potentialMap(pot)

      const MCF& const_mcf = mcf;

      const typename MCF::FlowMap &fm = const_mcf.flowMap();

      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();

      v = const_mcf.totalCost();

      double x = const_mcf.template totalCost<double>();

      v = const_mcf.flow(a);

      v = const_mcf.potential(n);

                const SM& supply, const FM& flow,

                ProblemType type = EQ )

    bool b = (type ==  EQ && sum == supply[n]) ||

             (type == GEQ && sum >= supply[n]) ||

             (type == LEQ && sum <= supply[n]);

    if (!b) return false;

           typename CM, typename SM, typename FM, typename PM >

                     const CM& cost, const SM& supply, const FM& flow, 

                     const PM& pi )

  for (NodeIt n(gr); opt && n != INVALID; ++n) {

    typename SM::Value sum = 0;

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

      sum += flow[e];

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

      sum -= flow[e];

    opt = (sum == supply[n]) || (pi[n] == 0);

  }

               const std::string &test_id = "",

               ProblemType type = EQ )

    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),

    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),

  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);

    .nodeMap("sup4", s4)

    .nodeMap("sup5", s5)

    // Check the equality form

    // Check the GEQ form

    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);

    checkMcf(mcf, mcf.run(),

             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);

    mcf.problemType(mcf.GEQ);

    checkMcf(mcf, mcf.lowerMap(l2).run(),

             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);

    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);

    checkMcf(mcf, mcf.run(),

             gr, l2, u, c, s5, false,    0, "#A11", GEQ);

    // Check the LEQ form

    mcf.reset().problemType(mcf.LEQ);

    mcf.upperMap(u).costMap(c).supplyMap(s5);

    checkMcf(mcf, mcf.run(),

             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);

    checkMcf(mcf, mcf.lowerMap(l2).run(),

             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);

    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);

    checkMcf(mcf, mcf.run(),

             gr, l2, u, c, s4, false,    0, "#A14", LEQ);