RhodeCode

Formally, let \f$G=(V,A)\f$ be a digraph, \f$lower: A\rightarrow\mathbf{Z}\f$,

\f$upper: A\rightarrow\mathbf{Z}\cup\{+\infty\}\f$ denote the lower and

\f$lower(uv) \leq upper(uv)\f$ must hold for all \f$uv\in A\f$,

\f$cost: A\rightarrow\mathbf{Z}\f$ denotes the cost per unit flow

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

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

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

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

\f$uv\in A\f$ with respect to the potential function \f$\pi\f$, i.e.

All algorithms provide dual solution (node potentials) as well,

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

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

  ///

  /// Most of the parameters of the problem (except for the digraph)

  /// can be given using separate functions, and the algorithm can be

  /// executed using the \ref run() function. If some parameters are not

  /// specified, then default values will be used.

    /// \brief Problem type constants for the \c run() function.

    /// Enum type containing the problem type constants that can be

    /// returned by the \ref run() function of the algorithm.

    enum ProblemType {

      /// The problem has no feasible solution (flow).

      INFEASIBLE,

      /// The problem has optimal solution (i.e. it is feasible and

      /// bounded), and the algorithm has found optimal flow and node

      /// potentials (primal and dual solutions).

      OPTIMAL,

      /// The objective function of the problem is unbounded, i.e.

      /// there is a directed cycle having negative total cost and

      /// infinite upper bound.

      UNBOUNDED

    };

    /// \brief Constants for selecting the type of the supply constraints.

    ///

    /// Enum type containing constants for selecting the supply 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 supply 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 supplyType()

    /// Note that the equality form is a special case of both supply types.

    enum SupplyType {

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

      /// supply/demand constraints in the definition, 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>

      /// supply/demand constraints in the definition, 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

    };

    /// \brief Constants for selecting the pivot rule.

    ///

    /// Enum type containing constants for selecting the pivot rule for

    /// the \ref run() function.

    ///

    SupplyType _stype;

    Flow _sum_supply;

  public:

    /// \brief Constant for infinite upper bounds (capacities).

    ///

    /// Constant for infinite upper bounds (capacities).

    /// It is \c std::numeric_limits<Flow>::infinity() if available,

    /// \c std::numeric_limits<Flow>::max() otherwise.

    const Flow INF;

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

      _node_id(graph),

      INF(std::numeric_limits<Flow>::has_infinity ?

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

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

      // Check the value types

      LEMON_ASSERT(std::numeric_limits<Flow>::is_signed,

        "The flow type of NetworkSimplex must be signed");

      LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,

        "The cost type of NetworkSimplex must be signed");

    /// If it is not used before calling \ref run(), the lower bounds

    /// will be set to zero on all arcs.

    template <typename LowerMap>

    NetworkSimplex& lowerMap(const LowerMap& map) {

    /// If it is not used before calling \ref run(), the upper bounds

    /// will be set to \ref INF on all arcs (i.e. the flow value will be

    /// unbounded from above on each arc).

    template<typename UpperMap>

    NetworkSimplex& upperMap(const UpperMap& map) {

    template<typename CostMap>

    NetworkSimplex& costMap(const CostMap& map) {

    template<typename SupplyMap>

    NetworkSimplex& supplyMap(const SupplyMap& map) {

    /// Using this function has the same effect as using \ref supplyMap()

    /// with such a map in which \c k is assigned to \c s, \c -k is

    /// assigned to \c t and all other nodes have zero supply value.

    ///

    /// \brief Set the type of the supply constraints.

    /// This function sets the type of the supply/demand constraints.

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

    /// For more information see \ref SupplyType.

    NetworkSimplex& supplyType(SupplyType supply_type) {

      _stype = supply_type;

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

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

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

    /// \return \c INFEASIBLE if no feasible flow exists,

    /// \n \c OPTIMAL if the problem has optimal solution

    /// (i.e. it is feasible and bounded), and the algorithm has found

    /// optimal flow and node potentials (primal and dual solutions),

    /// \n \c UNBOUNDED if the objective function of the problem is

    /// unbounded, i.e. there is a directed cycle having negative total

    /// cost and infinite upper bound.

    ///

    /// \see ProblemType, PivotRule

    ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) {

      if (!init()) return INFEASIBLE;

      return start(pivot_rule);

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

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

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

      _stype = GEQ;

    /// Its complexity is O(e).

    template <typename Value>

    Value totalCost() const {

      Value c = 0;

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

      _sum_supply = 0;

          _sum_supply += _supply[i];

        valid_supply = (_stype == GEQ && _sum_supply <= 0) ||

                       (_stype == LEQ && _sum_supply >= 0);

      Cost ART_COST;

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

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

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

        ART_COST = (ART_COST + 1) * _node_num;

      _supply[_root] = -_sum_supply;

      if (_sum_supply < 0) {

        _pi[_root] = -ART_COST;

        _pi[_root] = ART_COST;

            _cap[i] = INF;

          if (c > 0) {

            if (_cap[i] < INF) _cap[i] -= c;

            _supply[_source[i]] -= c;

            _supply[_target[i]] += c;

          }

          else if (c < 0) {

            if (_cap[i] < INF + c) {

              _cap[i] -= c;

            } else {

              _cap[i] = INF;

            }

        _cost[e] = ART_COST;

        _cap[e] = INF;

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

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

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

        d = _forward[u] ?

          _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]);

        d = _forward[u] ? 

          (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e];

    ProblemType start(PivotRule pivot_rule) {

      return INFEASIBLE; // avoid warning

    ProblemType start() {

        if (delta >= INF) return UNBOUNDED;

      // Check feasibility

      if (_sum_supply < 0) {

        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

          if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;

        }

      }

      else if (_sum_supply > 0) {

        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

          if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;

        }

      }

      else {

        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {

          if (_flow[e] != 0) return INFEASIBLE;

        }

      }

      return OPTIMAL;

#include <limits>

  "label  sup1 sup2 sup3 sup4 sup5 sup6\n"

  "    1    20   27    0   30   20   30\n"

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

  "    3     0    0    0    0    0    0\n"

  "    4     0    0    0    0    0    0\n"

  "    5     9    0    0    0    6   11\n"

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

  "    7     0    0    0    0    0    0\n"

  "    8     0    0    0    0    0    3\n"

  "    9     3    0    0    0    0    0\n"

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

  "   11     0    0    0    0  -10    0\n"

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

  "\n"                

  "       cost  cap low1 low2 low3\n"

  " 1  2    70   11    0    8    8\n"

  " 1  3   150    3    0    1    0\n"

  " 1  4    80   15    0    2    2\n"

  " 2  8    80   12    0    0    0\n"

  " 3  5   140    5    0    3    1\n"

  " 4  6    60   10    0    1    0\n"

  " 4  7    80    2    0    0    0\n"

  " 4  8   110    3    0    0    0\n"

  " 5  7    60   14    0    0    0\n"

  " 5 11   120   12    0    0    0\n"

  " 6  3     0    3    0    0    0\n"

  " 6  9   140    4    0    0    0\n"

  " 6 10    90    8    0    0    0\n"

  " 7  1    30    5    0    0   -5\n"

  " 8 12    60   16    0    4    3\n"

  " 9 12    50    6    0    0    0\n"

  "10 12    70   13    0    5    2\n"

  "10  2   100    7    0    0    0\n"

  "10  7    60   10    0    0   -3\n"

  "11 10    20   14    0    6  -20\n"

  "12 11    30   10    0    0  -10\n"

enum SupplyType {

      c = const_mcf.totalCost();

      c = const_mcf.potential(n);

      v = const_mcf.INF;

    Cost c;

                SupplyType type = EQ )

           typename CM, typename SM,

           typename PT >

void checkMcf( const MCF& mcf, PT mcf_result,

               PT result, bool optimal, typename CM::Value total,

               SupplyType type = EQ )

  if (optimal) {

  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), l3(gr), u(gr);

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

    .arcMap("low3", l3)

    .nodeMap("sup6", s6)

  // Build a test digraph for testing negative costs

  Digraph ngr;

  Node n1 = ngr.addNode();

  Node n2 = ngr.addNode();

  Node n3 = ngr.addNode();

  Node n4 = ngr.addNode();

  Node n5 = ngr.addNode();

  Node n6 = ngr.addNode();

  Node n7 = ngr.addNode();

  Arc a1 = ngr.addArc(n1, n2);

  Arc a2 = ngr.addArc(n1, n3);

  Arc a3 = ngr.addArc(n2, n4);

  Arc a4 = ngr.addArc(n3, n4);

  Arc a5 = ngr.addArc(n3, n2);

  Arc a6 = ngr.addArc(n5, n3);

  Arc a7 = ngr.addArc(n5, n6);

  Arc a8 = ngr.addArc(n6, n7);

  Arc a9 = ngr.addArc(n7, n5);

  Digraph::ArcMap<int> nc(ngr), nl1(ngr, 0), nl2(ngr, 0);

  ConstMap<Arc, int> nu1(std::numeric_limits<int>::max()), nu2(5000);

  Digraph::NodeMap<int> ns(ngr, 0);

  nl2[a7] =  1000;

  nl2[a8] = -1000;

  ns[n1] =  100;

  ns[n4] = -100;

  nc[a1] =  100;

  nc[a2] =   30;

  nc[a3] =   20;

  nc[a4] =   80;

  nc[a5] =   50;

  nc[a6] =   10;

  nc[a7] =   80;

  nc[a8] =   30;

  nc[a9] = -120;

             gr, l1, u, c, s1, mcf.OPTIMAL, true,   5240, "#A1");

             gr, l1, u, c, s2, mcf.OPTIMAL, true,   7620, "#A2");

             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#A3");

             gr, l2, u, c, s2, mcf.OPTIMAL, true,   8010, "#A4");

             gr, l1, cu, cc, s1, mcf.OPTIMAL, true,   74, "#A5");

             gr, l2, cu, cc, s2, mcf.OPTIMAL, true,   94, "#A6");

             gr, l1, cu, cc, s3, mcf.OPTIMAL, true,    0, "#A7");

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

             gr, l2, u, cc, s3, mcf.INFEASIBLE, false, 0, "#A8");

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

    checkMcf(mcf, mcf.run(),

             gr, l3, u, c, s4, mcf.OPTIMAL, true,   6360, "#A9");

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

             gr, l1, u, c, s5, mcf.OPTIMAL, true,   3530, "#A10", GEQ);

    mcf.supplyType(mcf.GEQ);

             gr, l2, u, c, s5, mcf.OPTIMAL, true,   4540, "#A11", GEQ);

    mcf.supplyType(mcf.CARRY_SUPPLIES).supplyMap(s6);

             gr, l2, u, c, s6, mcf.INFEASIBLE, false,  0, "#A12", GEQ);

    mcf.reset().supplyType(mcf.LEQ);

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

             gr, l1, u, c, s6, mcf.OPTIMAL, true,   5080, "#A13", LEQ);

             gr, l2, u, c, s6, mcf.OPTIMAL, true,   5930, "#A14", LEQ);

    mcf.supplyType(mcf.SATISFY_DEMANDS).supplyMap(s5);

             gr, l2, u, c, s5, mcf.INFEASIBLE, false,  0, "#A15", LEQ);

    // Check negative costs

    NetworkSimplex<Digraph> nmcf(ngr);

    nmcf.lowerMap(nl1).costMap(nc).supplyMap(ns);

    checkMcf(nmcf, nmcf.run(),

      ngr, nl1, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A16");

    checkMcf(nmcf, nmcf.upperMap(nu2).run(),

      ngr, nl1, nu2, nc, ns, nmcf.OPTIMAL, true,  -40000, "#A17");

    nmcf.reset().lowerMap(nl2).costMap(nc).supplyMap(ns);

    checkMcf(nmcf, nmcf.run(),

      ngr, nl2, nu1, nc, ns, nmcf.UNBOUNDED, false,    0, "#A18");

    mcf.supplyMap(s1).costMap(c).upperMap(u).lowerMap(l2);

             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B1");

             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B2");

             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B3");

             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B4");

             gr, l2, u, c, s1, mcf.OPTIMAL, true,   5970, "#B5");

  ns.lowerMap(lower).upperMap(cap).costMap(cost).supplyMap(sup);

  if (sum_sup > 0) ns.supplyType(ns.LEQ);