COIN-OR::LEMON - Graph Library

Changeset 656:e6927fe719e6 in lemon for lemon/network_simplex.h


Ignore:
Timestamp:
04/17/09 18:04:36 (10 years ago)
Author:
Peter Kovacs <kpeter@…>
Branch:
default
Phase:
public
Message:

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.

File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/network_simplex.h

    r655 r656  
    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 {
     
    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.
     
    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
     
    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:
     
    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:
     
    156221    Node _psource, _ptarget;
    157222    Flow _pstflow;
     223    ProblemType _ptype;
    158224
    159225    // Result maps
     
    587653    /// \brief Constructor.
    588654    ///
    589     /// Constructor.
     655    /// The constructor of the class.
    590656    ///
    591657    /// \param graph The digraph the algorithm runs on.
     
    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),
     
    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    ///
     
    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.
     
    796882      return *this;
    797883    }
     884   
     885    /// @}
    798886
    799887    /// \name Execution Control
     
    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);
     
    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
     
    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    }
     
    10011099      // Initialize node related data
    10021100      bool valid_supply = true;
     1101      Flow sum_supply = 0;
    10031102      if (!_pstsup && !_psupply) {
    10041103        _pstsup = true;
     
    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];
    1015         }
    1016         valid_supply = (sum == 0);
     1112          sum_supply += _supply[i];
     1113        }
     1114        valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
     1115                       (_ptype == LEQ && sum_supply >= 0);
    10171116      } else {
    10181117        int 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;
     1126
     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;
    10271211
    10281212      // Set data for the artificial root node
     
    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
     
    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) {
     
    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) {
     
    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      }
     
    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) {
Note: See TracChangeset for help on using the changeset viewer.