COIN-OR::LEMON - Graph Library

Ticket #219: ns-geq-leq-e6927fe719e6.patch

File ns-geq-leq-e6927fe719e6.patch, 28.1 KB (added by Peter Kovacs, 10 years ago)
  • doc/groups.dox

    # HG changeset patch
    # User Peter Kovacs <kpeter@inf.elte.hu>
    # Date 1239984276 -7200
    # Node ID e6927fe719e6e2cc97345c5f89f1cb0d66d2d0c9
    # Parent  6ac5d9ae1d3dfdfcc48802e03bc2160d5881a78c
    Support >= and <= constraints in NetworkSimplex (#219, #234)
    
    By default the same inequality constraints are supported as by
    Circulation (the GEQ form), but the LEQ form can also be selected
    using the problemType() function.
    
    The documentation of the min. cost flow module is reworked and
    extended with important notes and explanations about the different
    variants of the problem and about the dual solution and optimality
    conditions.
    
    diff --git a/doc/groups.dox b/doc/groups.dox
    a b  
    317317
    318318The \e maximum \e flow \e problem is to find a flow of maximum value between
    319319a single source and a single target. Formally, there is a \f$G=(V,A)\f$
    320 digraph, a \f$cap:A\rightarrow\mathbf{R}^+_0\f$ capacity function and
     320digraph, a \f$cap: A\rightarrow\mathbf{R}^+_0\f$ capacity function and
    321321\f$s, t \in V\f$ source and target nodes.
    322 A maximum flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of the
     322A maximum flow is an \f$f: A\rightarrow\mathbf{R}^+_0\f$ solution of the
    323323following optimization problem.
    324324
    325 \f[ \max\sum_{a\in\delta_{out}(s)}f(a) - \sum_{a\in\delta_{in}(s)}f(a) \f]
    326 \f[ \sum_{a\in\delta_{out}(v)} f(a) = \sum_{a\in\delta_{in}(v)} f(a)
    327     \qquad \forall v\in V\setminus\{s,t\} \f]
    328 \f[ 0 \leq f(a) \leq cap(a) \qquad \forall a\in A \f]
     325\f[ \max\sum_{sv\in A} f(sv) - \sum_{vs\in A} f(vs) \f]
     326\f[ \sum_{uv\in A} f(uv) = \sum_{vu\in A} f(vu)
     327    \quad \forall u\in V\setminus\{s,t\} \f]
     328\f[ 0 \leq f(uv) \leq cap(uv) \quad \forall uv\in A \f]
    329329
    330330LEMON contains several algorithms for solving maximum flow problems:
    331331- \ref EdmondsKarp Edmonds-Karp algorithm.
     
    345345
    346346\brief Algorithms for finding minimum cost flows and circulations.
    347347
    348 This group describes the algorithms for finding minimum cost flows and
     348This group contains the algorithms for finding minimum cost flows and
    349349circulations.
    350350
    351351The \e minimum \e cost \e flow \e problem is to find a feasible flow of
    352352minimum total cost from a set of supply nodes to a set of demand nodes
    353 in a network with capacity constraints and arc costs.
     353in a network with capacity constraints (lower and upper bounds)
     354and arc costs.
    354355Formally, let \f$G=(V,A)\f$ be a digraph,
    355356\f$lower, upper: A\rightarrow\mathbf{Z}^+_0\f$ denote the lower and
    356 upper bounds for the flow values on the arcs,
     357upper bounds for the flow values on the arcs, for which
     358\f$0 \leq lower(uv) \leq upper(uv)\f$ holds for all \f$uv\in A\f$.
    357359\f$cost: A\rightarrow\mathbf{Z}^+_0\f$ denotes the cost per unit flow
    358 on the arcs, and
    359 \f$supply: V\rightarrow\mathbf{Z}\f$ denotes the supply/demand values
    360 of the nodes.
    361 A minimum cost flow is an \f$f:A\rightarrow\mathbf{R}^+_0\f$ solution of
    362 the following optimization problem.
     360on the arcs, and \f$sup: V\rightarrow\mathbf{Z}\f$ denotes the
     361signed supply values of the nodes.
     362If \f$sup(u)>0\f$, then \f$u\f$ is a supply node with \f$sup(u)\f$
     363supply, if \f$sup(u)<0\f$, then \f$u\f$ is a demand node with
     364\f$-sup(u)\f$ demand.
     365A minimum cost flow is an \f$f: A\rightarrow\mathbf{Z}^+_0\f$ solution
     366of the following optimization problem.
    363367
    364 \f[ \min\sum_{a\in A} f(a) cost(a) \f]
    365 \f[ \sum_{a\in\delta_{out}(v)} f(a) - \sum_{a\in\delta_{in}(v)} f(a) =
    366     supply(v) \qquad \forall v\in V \f]
    367 \f[ lower(a) \leq f(a) \leq upper(a) \qquad \forall a\in A \f]
     368\f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
     369\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
     370    sup(u) \quad \forall u\in V \f]
     371\f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    368372
    369 LEMON contains several algorithms for solving minimum cost flow problems:
    370  - \ref CycleCanceling Cycle-canceling algorithms.
    371  - \ref CapacityScaling Successive shortest path algorithm with optional
     373The sum of the supply values, i.e. \f$\sum_{u\in V} sup(u)\f$ must be
     374zero or negative in order to have a feasible solution (since the sum
     375of the expressions on the left-hand side of the inequalities is zero).
     376It means that the total demand must be greater or equal to the total
     377supply and all the supplies have to be carried out from the supply nodes,
     378but there could be demands that are not satisfied.
     379If \f$\sum_{u\in V} sup(u)\f$ is zero, then all the supply/demand
     380constraints have to be satisfied with equality, i.e. all demands
     381have to be satisfied and all supplies have to be used.
     382
     383If you need the opposite inequalities in the supply/demand constraints
     384(i.e. the total demand is less than the total supply and all the demands
     385have to be satisfied while there could be supplies that are not used),
     386then you could easily transform the problem to the above form by reversing
     387the direction of the arcs and taking the negative of the supply values
     388(e.g. using \ref ReverseDigraph and \ref NegMap adaptors).
     389However \ref NetworkSimplex algorithm also supports this form directly
     390for the sake of convenience.
     391
     392A feasible solution for this problem can be found using \ref Circulation.
     393
     394Note that the above formulation is actually more general than the usual
     395definition of the minimum cost flow problem, in which strict equalities
     396are required in the supply/demand contraints, i.e.
     397
     398\f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) =
     399    sup(u) \quad \forall u\in V. \f]
     400
     401However if the sum of the supply values is zero, then these two problems
     402are equivalent. So if you need the equality form, you have to ensure this
     403additional contraint for the algorithms.
     404
     405The dual solution of the minimum cost flow problem is represented by node
     406potentials \f$\pi: V\rightarrow\mathbf{Z}\f$.
     407An \f$f: A\rightarrow\mathbf{Z}^+_0\f$ feasible solution of the problem
     408is optimal if and only if for some \f$\pi: V\rightarrow\mathbf{Z}\f$
     409node potentials the following \e complementary \e slackness optimality
     410conditions hold.
     411
     412 - For all \f$uv\in A\f$ arcs:
     413   - if \f$cost^\pi(uv)>0\f$, then \f$f(uv)=lower(uv)\f$;
     414   - if \f$lower(uv)<f(uv)<upper(uv)\f$, then \f$cost^\pi(uv)=0\f$;
     415   - if \f$cost^\pi(uv)<0\f$, then \f$f(uv)=upper(uv)\f$.
     416 - For all \f$u\in V\f$:
     417   - if \f$\sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \neq sup(u)\f$,
     418     then \f$\pi(u)=0\f$.
     419 
     420Here \f$cost^\pi(uv)\f$ denotes the \e reduced \e cost of the arc
     421\f$uv\in A\f$ with respect to the node potentials \f$\pi\f$, i.e.
     422\f[ cost^\pi(uv) = cost(uv) + \pi(u) - \pi(v).\f]
     423
     424All algorithms provide dual solution (node potentials) as well
     425if an optimal flow is found.
     426
     427LEMON contains several algorithms for solving minimum cost flow problems.
     428 - \ref NetworkSimplex Primal Network Simplex algorithm with various
     429   pivot strategies.
     430 - \ref CostScaling Push-Relabel and Augment-Relabel algorithms based on
     431   cost scaling.
     432 - \ref CapacityScaling Successive Shortest %Path algorithm with optional
    372433   capacity scaling.
    373  - \ref CostScaling Push-relabel and augment-relabel algorithms based on
    374    cost scaling.
    375  - \ref NetworkSimplex Primal network simplex algorithm with various
    376    pivot strategies.
     434 - \ref CancelAndTighten The Cancel and Tighten algorithm.
     435 - \ref CycleCanceling Cycle-Canceling algorithms.
     436
     437Most of these implementations support the general inequality form of the
     438minimum cost flow problem, but CancelAndTighten and CycleCanceling
     439only support the equality form due to the primal method they use.
     440
     441In general NetworkSimplex is the most efficient implementation,
     442but in special cases other algorithms could be faster.
     443For example, if the total supply and/or capacities are rather small,
     444CapacityScaling is usually the fastest algorithm (without effective scaling).
    377445*/
    378446
    379447/**
  • lemon/network_simplex.h

    diff --git a/lemon/network_simplex.h b/lemon/network_simplex.h
    a b  
    3030
    3131#include <lemon/core.h>
    3232#include <lemon/math.h>
     33#include <lemon/maps.h>
     34#include <lemon/circulation.h>
     35#include <lemon/adaptors.h>
    3336
    3437namespace lemon {
    3538
     
    4750  ///
    4851  /// In general this class is the fastest implementation available
    4952  /// in LEMON for the minimum cost flow problem.
     53  /// Moreover it supports both direction of the supply/demand inequality
     54  /// constraints. For more information see \ref ProblemType.
    5055  ///
    5156  /// \tparam GR The digraph type the algorithm runs on.
    5257  /// \tparam F The value type used for flow amounts, capacity bounds
     
    5863  /// be integer.
    5964  ///
    6065  /// \note %NetworkSimplex provides five different pivot rule
    61   /// implementations. For more information see \ref PivotRule.
     66  /// implementations, from which the most efficient one is used
     67  /// by default. For more information see \ref PivotRule.
    6268  template <typename GR, typename F = int, typename C = F>
    6369  class NetworkSimplex
    6470  {
     
    6874    typedef F Flow;
    6975    /// The cost type of the algorithm
    7076    typedef C Cost;
     77#ifdef DOXYGEN
     78    /// The type of the flow map
     79    typedef GR::ArcMap<Flow> FlowMap;
     80    /// The type of the potential map
     81    typedef GR::NodeMap<Cost> PotentialMap;
     82#else
    7183    /// The type of the flow map
    7284    typedef typename GR::template ArcMap<Flow> FlowMap;
    7385    /// The type of the potential map
    7486    typedef typename GR::template NodeMap<Cost> PotentialMap;
     87#endif
    7588
    7689  public:
    7790
     
    117130      /// candidate list and extends this list in every iteration.
    118131      ALTERING_LIST
    119132    };
     133   
     134    /// \brief Enum type for selecting the problem type.
     135    ///
     136    /// Enum type for selecting the problem type, i.e. the direction of
     137    /// the inequalities in the supply/demand constraints of the
     138    /// \ref min_cost_flow "minimum cost flow problem".
     139    ///
     140    /// The default problem type is \c GEQ, since this form is supported
     141    /// by other minimum cost flow algorithms and the \ref Circulation
     142    /// algorithm as well.
     143    /// The \c LEQ problem type can be selected using the \ref problemType()
     144    /// function.
     145    ///
     146    /// Note that the equality form is a special case of both problem type.
     147    enum ProblemType {
     148
     149      /// This option means that there are "<em>greater or equal</em>"
     150      /// constraints in the defintion, i.e. the exact formulation of the
     151      /// problem is the following.
     152      /**
     153          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
     154          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
     155              sup(u) \quad \forall u\in V \f]
     156          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
     157      */
     158      /// It means that the total demand must be greater or equal to the
     159      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
     160      /// negative) and all the supplies have to be carried out from
     161      /// the supply nodes, but there could be demands that are not
     162      /// satisfied.
     163      GEQ,
     164      /// It is just an alias for the \c GEQ option.
     165      CARRY_SUPPLIES = GEQ,
     166
     167      /// This option means that there are "<em>less or equal</em>"
     168      /// constraints in the defintion, i.e. the exact formulation of the
     169      /// problem is the following.
     170      /**
     171          \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
     172          \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
     173              sup(u) \quad \forall u\in V \f]
     174          \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
     175      */
     176      /// It means that the total demand must be less or equal to the
     177      /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
     178      /// positive) and all the demands have to be satisfied, but there
     179      /// could be supplies that are not carried out from the supply
     180      /// nodes.
     181      LEQ,
     182      /// It is just an alias for the \c LEQ option.
     183      SATISFY_DEMANDS = LEQ
     184    };
    120185
    121186  private:
    122187
     
    155220    bool _pstsup;
    156221    Node _psource, _ptarget;
    157222    Flow _pstflow;
     223    ProblemType _ptype;
    158224
    159225    // Result maps
    160226    FlowMap *_flow_map;
     
    586652
    587653    /// \brief Constructor.
    588654    ///
    589     /// Constructor.
     655    /// The constructor of the class.
    590656    ///
    591657    /// \param graph The digraph the algorithm runs on.
    592658    NetworkSimplex(const GR& graph) :
    593659      _graph(graph),
    594660      _plower(NULL), _pupper(NULL), _pcost(NULL),
    595       _psupply(NULL), _pstsup(false),
     661      _psupply(NULL), _pstsup(false), _ptype(GEQ),
    596662      _flow_map(NULL), _potential_map(NULL),
    597663      _local_flow(false), _local_potential(false),
    598664      _node_id(graph)
     
    611677      if (_local_potential) delete _potential_map;
    612678    }
    613679
     680    /// \name Parameters
     681    /// The parameters of the algorithm can be specified using these
     682    /// functions.
     683
     684    /// @{
     685
    614686    /// \brief Set the lower bounds on the arcs.
    615687    ///
    616688    /// This function sets the lower bounds on the arcs.
     
    760832      _pstflow = k;
    761833      return *this;
    762834    }
     835   
     836    /// \brief Set the problem type.
     837    ///
     838    /// This function sets the problem type for the algorithm.
     839    /// If it is not used before calling \ref run(), the \ref GEQ problem
     840    /// type will be used.
     841    ///
     842    /// For more information see \ref ProblemType.
     843    ///
     844    /// \return <tt>(*this)</tt>
     845    NetworkSimplex& problemType(ProblemType problem_type) {
     846      _ptype = problem_type;
     847      return *this;
     848    }
    763849
    764850    /// \brief Set the flow map.
    765851    ///
     
    795881      _potential_map = &map;
    796882      return *this;
    797883    }
     884   
     885    /// @}
    798886
    799887    /// \name Execution Control
    800888    /// The algorithm can be executed using \ref run().
     
    804892    /// \brief Run the algorithm.
    805893    ///
    806894    /// This function runs the algorithm.
    807     /// The paramters can be specified using \ref lowerMap(),
     895    /// The paramters can be specified using functions \ref lowerMap(),
    808896    /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
    809     /// \ref costMap(), \ref supplyMap() and \ref stSupply()
    810     /// functions. For example,
     897    /// \ref costMap(), \ref supplyMap(), \ref stSupply(),
     898    /// \ref problemType(), \ref flowMap() and \ref potentialMap().
     899    /// For example,
    811900    /// \code
    812901    ///   NetworkSimplex<ListDigraph> ns(graph);
    813902    ///   ns.boundMaps(lower, upper).costMap(cost)
     
    830919    /// \brief Reset all the parameters that have been given before.
    831920    ///
    832921    /// This function resets all the paramaters that have been given
    833     /// using \ref lowerMap(), \ref upperMap(), \ref capacityMap(),
    834     /// \ref boundMaps(), \ref costMap(), \ref supplyMap() and
    835     /// \ref stSupply() functions before.
     922    /// before using functions \ref lowerMap(), \ref upperMap(),
     923    /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
     924    /// \ref supplyMap(), \ref stSupply(), \ref problemType(),
     925    /// \ref flowMap() and \ref potentialMap().
    836926    ///
    837927    /// It is useful for multiple run() calls. If this function is not
    838928    /// used, all the parameters given before are kept for the next
     
    869959      _pcost = NULL;
    870960      _psupply = NULL;
    871961      _pstsup = false;
     962      _ptype = GEQ;
     963      if (_local_flow) delete _flow_map;
     964      if (_local_potential) delete _potential_map;
     965      _flow_map = NULL;
     966      _potential_map = NULL;
     967      _local_flow = false;
     968      _local_potential = false;
     969
    872970      return *this;
    873971    }
    874972
     
    10001098
    10011099      // Initialize node related data
    10021100      bool valid_supply = true;
     1101      Flow sum_supply = 0;
    10031102      if (!_pstsup && !_psupply) {
    10041103        _pstsup = true;
    10051104        _psource = _ptarget = NodeIt(_graph);
    10061105        _pstflow = 0;
    10071106      }
    10081107      if (_psupply) {
    1009         Flow sum = 0;
    10101108        int i = 0;
    10111109        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
    10121110          _node_id[n] = i;
    10131111          _supply[i] = (*_psupply)[n];
    1014           sum += _supply[i];
     1112          sum_supply += _supply[i];
    10151113        }
    1016         valid_supply = (sum == 0);
     1114        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
     1115                       (_ptype == LEQ && sum_supply >= 0);
    10171116      } else {
    10181117        int i = 0;
    10191118        for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
     
    10211120          _supply[i] = 0;
    10221121        }
    10231122        _supply[_node_id[_psource]] =  _pstflow;
    1024         _supply[_node_id[_ptarget]]   = -_pstflow;
     1123        _supply[_node_id[_ptarget]] = -_pstflow;
    10251124      }
    10261125      if (!valid_supply) return false;
    10271126
     1127      // Infinite capacity value
     1128      Flow inf_cap =
     1129        std::numeric_limits<Flow>::has_infinity ?
     1130        std::numeric_limits<Flow>::infinity() :
     1131        std::numeric_limits<Flow>::max();
     1132
     1133      // Initialize artifical cost
     1134      Cost art_cost;
     1135      if (std::numeric_limits<Cost>::is_exact) {
     1136        art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
     1137      } else {
     1138        art_cost = std::numeric_limits<Cost>::min();
     1139        for (int i = 0; i != _arc_num; ++i) {
     1140          if (_cost[i] > art_cost) art_cost = _cost[i];
     1141        }
     1142        art_cost = (art_cost + 1) * _node_num;
     1143      }
     1144
     1145      // Run Circulation to check if a feasible solution exists
     1146      typedef ConstMap<Arc, Flow> ConstArcMap;
     1147      FlowNodeMap *csup = NULL;
     1148      bool local_csup = false;
     1149      if (_psupply) {
     1150        csup = _psupply;
     1151      } else {
     1152        csup = new FlowNodeMap(_graph, 0);
     1153        (*csup)[_psource] =  _pstflow;
     1154        (*csup)[_ptarget] = -_pstflow;
     1155        local_csup = true;
     1156      }
     1157      bool circ_result = false;
     1158      if (_ptype == GEQ || (_ptype == LEQ && sum_supply == 0)) {
     1159        // GEQ problem type
     1160        if (_plower) {
     1161          if (_pupper) {
     1162            Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap>
     1163              circ(_graph, *_plower, *_pupper, *csup);
     1164            circ_result = circ.run();
     1165          } else {
     1166            Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap>
     1167              circ(_graph, *_plower, ConstArcMap(inf_cap), *csup);
     1168            circ_result = circ.run();
     1169          }
     1170        } else {
     1171          if (_pupper) {
     1172            Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap>
     1173              circ(_graph, ConstArcMap(0), *_pupper, *csup);
     1174            circ_result = circ.run();
     1175          } else {
     1176            Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap>
     1177              circ(_graph, ConstArcMap(0), ConstArcMap(inf_cap), *csup);
     1178            circ_result = circ.run();
     1179          }
     1180        }
     1181      } else {
     1182        // LEQ problem type
     1183        typedef ReverseDigraph<const GR> RevGraph;
     1184        typedef NegMap<FlowNodeMap> NegNodeMap;
     1185        RevGraph rgraph(_graph);
     1186        NegNodeMap neg_csup(*csup);
     1187        if (_plower) {
     1188          if (_pupper) {
     1189            Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap>
     1190              circ(rgraph, *_plower, *_pupper, neg_csup);
     1191            circ_result = circ.run();
     1192          } else {
     1193            Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap>
     1194              circ(rgraph, *_plower, ConstArcMap(inf_cap), neg_csup);
     1195            circ_result = circ.run();
     1196          }
     1197        } else {
     1198          if (_pupper) {
     1199            Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap>
     1200              circ(rgraph, ConstArcMap(0), *_pupper, neg_csup);
     1201            circ_result = circ.run();
     1202          } else {
     1203            Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap>
     1204              circ(rgraph, ConstArcMap(0), ConstArcMap(inf_cap), neg_csup);
     1205            circ_result = circ.run();
     1206          }
     1207        }
     1208      }
     1209      if (local_csup) delete csup;
     1210      if (!circ_result) return false;
     1211
    10281212      // Set data for the artificial root node
    10291213      _root = _node_num;
    10301214      _parent[_root] = -1;
     
    10331217      _rev_thread[0] = _root;
    10341218      _succ_num[_root] = all_node_num;
    10351219      _last_succ[_root] = _root - 1;
    1036       _supply[_root] = 0;
    1037       _pi[_root] = 0;
     1220      _supply[_root] = -sum_supply;
     1221      if (sum_supply < 0) {
     1222        _pi[_root] = -art_cost;
     1223      } else {
     1224        _pi[_root] = art_cost;
     1225      }
    10381226
    10391227      // Store the arcs in a mixed order
    10401228      int k = std::max(int(sqrt(_arc_num)), 10);
     
    10451233      }
    10461234
    10471235      // Initialize arc maps
    1048       Flow inf_cap =
    1049         std::numeric_limits<Flow>::has_infinity ?
    1050         std::numeric_limits<Flow>::infinity() :
    1051         std::numeric_limits<Flow>::max();
    10521236      if (_pupper && _pcost) {
    10531237        for (int i = 0; i != _arc_num; ++i) {
    10541238          Arc e = _arc_ref[i];
     
    10831267        }
    10841268      }
    10851269     
    1086       // Initialize artifical cost
    1087       Cost art_cost;
    1088       if (std::numeric_limits<Cost>::is_exact) {
    1089         art_cost = std::numeric_limits<Cost>::max() / 4 + 1;
    1090       } else {
    1091         art_cost = std::numeric_limits<Cost>::min();
    1092         for (int i = 0; i != _arc_num; ++i) {
    1093           if (_cost[i] > art_cost) art_cost = _cost[i];
    1094         }
    1095         art_cost = (art_cost + 1) * _node_num;
    1096       }
    1097 
    10981270      // Remove non-zero lower bounds
    10991271      if (_plower) {
    11001272        for (int i = 0; i != _arc_num; ++i) {
     
    11181290        _cost[e] = art_cost;
    11191291        _cap[e] = inf_cap;
    11201292        _state[e] = STATE_TREE;
    1121         if (_supply[u] >= 0) {
     1293        if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
    11221294          _flow[e] = _supply[u];
    11231295          _forward[u] = true;
    1124           _pi[u] = -art_cost;
     1296          _pi[u] = -art_cost + _pi[_root];
    11251297        } else {
    11261298          _flow[e] = -_supply[u];
    11271299          _forward[u] = false;
    1128           _pi[u] = art_cost;
     1300          _pi[u] = art_cost + _pi[_root];
    11291301        }
    11301302      }
    11311303
     
    13821554        }
    13831555      }
    13841556
    1385       // Check if the flow amount equals zero on all the artificial arcs
    1386       for (int e = _arc_num; e != _arc_num + _node_num; ++e) {
    1387         if (_flow[e] > 0) return false;
    1388       }
    1389 
    13901557      // Copy flow values to _flow_map
    13911558      if (_plower) {
    13921559        for (int i = 0; i != _arc_num; ++i) {
  • test/min_cost_flow_test.cc

    diff --git a/test/min_cost_flow_test.cc b/test/min_cost_flow_test.cc
    a b  
    3333
    3434char test_lgf[] =
    3535  "@nodes\n"
    36   "label  sup1 sup2 sup3\n"
    37   "    1    20   27    0\n"
    38   "    2    -4    0    0\n"
    39   "    3     0    0    0\n"
    40   "    4     0    0    0\n"
    41   "    5     9    0    0\n"
    42   "    6    -6    0    0\n"
    43   "    7     0    0    0\n"
    44   "    8     0    0    0\n"
    45   "    9     3    0    0\n"
    46   "   10    -2    0    0\n"
    47   "   11     0    0    0\n"
    48   "   12   -20  -27    0\n"
     36  "label  sup1 sup2 sup3 sup4 sup5\n"
     37  "    1    20   27    0   20   30\n"
     38  "    2    -4    0    0   -8   -3\n"
     39  "    3     0    0    0    0    0\n"
     40  "    4     0    0    0    0    0\n"
     41  "    5     9    0    0    6   11\n"
     42  "    6    -6    0    0   -5   -6\n"
     43  "    7     0    0    0    0    0\n"
     44  "    8     0    0    0    0    3\n"
     45  "    9     3    0    0    0    0\n"
     46  "   10    -2    0    0   -7   -2\n"
     47  "   11     0    0    0  -10    0\n"
     48  "   12   -20  -27    0  -30  -20\n"
    4949  "\n"
    5050  "@arcs\n"
    5151  "       cost  cap low1 low2\n"
     
    7676  "target 12\n";
    7777
    7878
     79enum ProblemType {
     80  EQ,
     81  GEQ,
     82  LEQ
     83};
     84
    7985// Check the interface of an MCF algorithm
    8086template <typename GR, typename Flow, typename Cost>
    8187class McfClassConcept
     
    97103             .costMap(cost)
    98104             .supplyMap(sup)
    99105             .stSupply(n, n, k)
     106             .flowMap(flow)
     107             .potentialMap(pot)
    100108             .run();
     109     
     110      const MCF& const_mcf = mcf;
    101111
    102       const typename MCF::FlowMap &fm = mcf.flowMap();
    103       const typename MCF::PotentialMap &pm = mcf.potentialMap();
     112      const typename MCF::FlowMap &fm = const_mcf.flowMap();
     113      const typename MCF::PotentialMap &pm = const_mcf.potentialMap();
    104114
    105       v = mcf.totalCost();
    106       double x = mcf.template totalCost<double>();
    107       v = mcf.flow(a);
    108       v = mcf.potential(n);
    109       mcf.flowMap(flow);
    110       mcf.potentialMap(pot);
     115      v = const_mcf.totalCost();
     116      double x = const_mcf.template totalCost<double>();
     117      v = const_mcf.flow(a);
     118      v = const_mcf.potential(n);
    111119
    112120      ignore_unused_variable_warning(fm);
    113121      ignore_unused_variable_warning(pm);
     
    142150template < typename GR, typename LM, typename UM,
    143151           typename SM, typename FM >
    144152bool checkFlow( const GR& gr, const LM& lower, const UM& upper,
    145                 const SM& supply, const FM& flow )
     153                const SM& supply, const FM& flow,
     154                ProblemType type = EQ )
    146155{
    147156  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    148157
     
    156165      sum += flow[e];
    157166    for (InArcIt e(gr, n); e != INVALID; ++e)
    158167      sum -= flow[e];
    159     if (sum != supply[n]) return false;
     168    bool b = (type ==  EQ && sum == supply[n]) ||
     169             (type == GEQ && sum >= supply[n]) ||
     170             (type == LEQ && sum <= supply[n]);
     171    if (!b) return false;
    160172  }
    161173
    162174  return true;
     
    165177// Check the feasibility of the given potentials (dual soluiton)
    166178// using the "Complementary Slackness" optimality condition
    167179template < typename GR, typename LM, typename UM,
    168            typename CM, typename FM, typename PM >
     180           typename CM, typename SM, typename FM, typename PM >
    169181bool checkPotential( const GR& gr, const LM& lower, const UM& upper,
    170                      const CM& cost, const FM& flow, const PM& pi )
     182                     const CM& cost, const SM& supply, const FM& flow,
     183                     const PM& pi )
    171184{
    172185  TEMPLATE_DIGRAPH_TYPEDEFS(GR);
    173186
     
    179192          (red_cost > 0 && flow[e] == lower[e]) ||
    180193          (red_cost < 0 && flow[e] == upper[e]);
    181194  }
     195 
     196  for (NodeIt n(gr); opt && n != INVALID; ++n) {
     197    typename SM::Value sum = 0;
     198    for (OutArcIt e(gr, n); e != INVALID; ++e)
     199      sum += flow[e];
     200    for (InArcIt e(gr, n); e != INVALID; ++e)
     201      sum -= flow[e];
     202    opt = (sum == supply[n]) || (pi[n] == 0);
     203  }
     204 
    182205  return opt;
    183206}
    184207
     
    190213               const GR& gr, const LM& lower, const UM& upper,
    191214               const CM& cost, const SM& supply,
    192215               bool result, typename CM::Value total,
    193                const std::string &test_id = "" )
     216               const std::string &test_id = "",
     217               ProblemType type = EQ )
    194218{
    195219  check(mcf_result == result, "Wrong result " + test_id);
    196220  if (result) {
    197     check(checkFlow(gr, lower, upper, supply, mcf.flowMap()),
     221    check(checkFlow(gr, lower, upper, supply, mcf.flowMap(), type),
    198222          "The flow is not feasible " + test_id);
    199223    check(mcf.totalCost() == total, "The flow is not optimal " + test_id);
    200     check(checkPotential(gr, lower, upper, cost, mcf.flowMap(),
     224    check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(),
    201225                         mcf.potentialMap()),
    202226          "Wrong potentials " + test_id);
    203227  }
     
    226250  // Read the test digraph
    227251  Digraph gr;
    228252  Digraph::ArcMap<int> c(gr), l1(gr), l2(gr), u(gr);
    229   Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr);
     253  Digraph::NodeMap<int> s1(gr), s2(gr), s3(gr), s4(gr), s5(gr);
    230254  ConstMap<Arc, int> cc(1), cu(std::numeric_limits<int>::max());
    231255  Node v, w;
    232256
     
    239263    .nodeMap("sup1", s1)
    240264    .nodeMap("sup2", s2)
    241265    .nodeMap("sup3", s3)
     266    .nodeMap("sup4", s4)
     267    .nodeMap("sup5", s5)
    242268    .node("source", v)
    243269    .node("target", w)
    244270    .run();
     
    247273  {
    248274    NetworkSimplex<Digraph> mcf(gr);
    249275
     276    // Check the equality form
    250277    mcf.upperMap(u).costMap(c);
    251278    checkMcf(mcf, mcf.supplyMap(s1).run(),
    252279             gr, l1, u, c, s1, true,  5240, "#A1");
     
    267294             gr, l1, cu, cc, s3, true,   0, "#A7");
    268295    checkMcf(mcf, mcf.boundMaps(l2, u).run(),
    269296             gr, l2, u, cc, s3, false,   0, "#A8");
     297
     298    // Check the GEQ form
     299    mcf.reset().upperMap(u).costMap(c).supplyMap(s4);
     300    checkMcf(mcf, mcf.run(),
     301             gr, l1, u, c, s4, true,  3530, "#A9", GEQ);
     302    mcf.problemType(mcf.GEQ);
     303    checkMcf(mcf, mcf.lowerMap(l2).run(),
     304             gr, l2, u, c, s4, true,  4540, "#A10", GEQ);
     305    mcf.problemType(mcf.CARRY_SUPPLIES).supplyMap(s5);
     306    checkMcf(mcf, mcf.run(),
     307             gr, l2, u, c, s5, false,    0, "#A11", GEQ);
     308
     309    // Check the LEQ form
     310    mcf.reset().problemType(mcf.LEQ);
     311    mcf.upperMap(u).costMap(c).supplyMap(s5);
     312    checkMcf(mcf, mcf.run(),
     313             gr, l1, u, c, s5, true,  5080, "#A12", LEQ);
     314    checkMcf(mcf, mcf.lowerMap(l2).run(),
     315             gr, l2, u, c, s5, true,  5930, "#A13", LEQ);
     316    mcf.problemType(mcf.SATISFY_DEMANDS).supplyMap(s4);
     317    checkMcf(mcf, mcf.run(),
     318             gr, l2, u, c, s4, false,    0, "#A14", LEQ);
    270319  }
    271320
    272321  // B. Test NetworkSimplex with each pivot rule