lemon/network_simplex.h
changeset 621 b536eaacb39b
parent 608 6ac5d9ae1d3d
child 612 0c8e5c688440
child 613 b1811c363299
equal deleted inserted replaced
6:9cdbdb665726 7:2a9d8e31081e
    28 #include <limits>
    28 #include <limits>
    29 #include <algorithm>
    29 #include <algorithm>
    30 
    30 
    31 #include <lemon/core.h>
    31 #include <lemon/core.h>
    32 #include <lemon/math.h>
    32 #include <lemon/math.h>
       
    33 #include <lemon/maps.h>
       
    34 #include <lemon/circulation.h>
       
    35 #include <lemon/adaptors.h>
    33 
    36 
    34 namespace lemon {
    37 namespace lemon {
    35 
    38 
    36   /// \addtogroup min_cost_flow
    39   /// \addtogroup min_cost_flow
    37   /// @{
    40   /// @{
    45   /// simplex method directly for the minimum cost flow problem.
    48   /// simplex method directly for the minimum cost flow problem.
    46   /// It is one of the most efficient solution methods.
    49   /// It is one of the most efficient solution methods.
    47   ///
    50   ///
    48   /// In general this class is the fastest implementation available
    51   /// In general this class is the fastest implementation available
    49   /// in LEMON for the minimum cost flow problem.
    52   /// 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.
    50   ///
    55   ///
    51   /// \tparam GR The digraph type the algorithm runs on.
    56   /// \tparam GR The digraph type the algorithm runs on.
    52   /// \tparam F The value type used for flow amounts, capacity bounds
    57   /// \tparam F The value type used for flow amounts, capacity bounds
    53   /// and supply values in the algorithm. By default it is \c int.
    58   /// and supply values in the algorithm. By default it is \c int.
    54   /// \tparam C The value type used for costs and potentials in the
    59   /// \tparam C The value type used for costs and potentials in the
    56   ///
    61   ///
    57   /// \warning Both value types must be signed and all input data must
    62   /// \warning Both value types must be signed and all input data must
    58   /// be integer.
    63   /// be integer.
    59   ///
    64   ///
    60   /// \note %NetworkSimplex provides five different pivot rule
    65   /// \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.
    62   template <typename GR, typename F = int, typename C = F>
    68   template <typename GR, typename F = int, typename C = F>
    63   class NetworkSimplex
    69   class NetworkSimplex
    64   {
    70   {
    65   public:
    71   public:
    66 
    72 
    67     /// The flow type of the algorithm
    73     /// The flow type of the algorithm
    68     typedef F Flow;
    74     typedef F Flow;
    69     /// The cost type of the algorithm
    75     /// The cost type of the algorithm
    70     typedef C Cost;
    76     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
    71     /// The type of the flow map
    83     /// The type of the flow map
    72     typedef typename GR::template ArcMap<Flow> FlowMap;
    84     typedef typename GR::template ArcMap<Flow> FlowMap;
    73     /// The type of the potential map
    85     /// The type of the potential map
    74     typedef typename GR::template NodeMap<Cost> PotentialMap;
    86     typedef typename GR::template NodeMap<Cost> PotentialMap;
       
    87 #endif
    75 
    88 
    76   public:
    89   public:
    77 
    90 
    78     /// \brief Enum type for selecting the pivot rule.
    91     /// \brief Enum type for selecting the pivot rule.
    79     ///
    92     ///
   115       /// It is a modified version of the Candidate List method.
   128       /// It is a modified version of the Candidate List method.
   116       /// It keeps only the several best eligible arcs from the former
   129       /// It keeps only the several best eligible arcs from the former
   117       /// candidate list and extends this list in every iteration.
   130       /// candidate list and extends this list in every iteration.
   118       ALTERING_LIST
   131       ALTERING_LIST
   119     };
   132     };
       
   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     };
   120 
   185 
   121   private:
   186   private:
   122 
   187 
   123     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   188     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   124 
   189 
   153     CostArcMap *_pcost;
   218     CostArcMap *_pcost;
   154     FlowNodeMap *_psupply;
   219     FlowNodeMap *_psupply;
   155     bool _pstsup;
   220     bool _pstsup;
   156     Node _psource, _ptarget;
   221     Node _psource, _ptarget;
   157     Flow _pstflow;
   222     Flow _pstflow;
       
   223     ProblemType _ptype;
   158 
   224 
   159     // Result maps
   225     // Result maps
   160     FlowMap *_flow_map;
   226     FlowMap *_flow_map;
   161     PotentialMap *_potential_map;
   227     PotentialMap *_potential_map;
   162     bool _local_flow;
   228     bool _local_flow;
   584 
   650 
   585   public:
   651   public:
   586 
   652 
   587     /// \brief Constructor.
   653     /// \brief Constructor.
   588     ///
   654     ///
   589     /// Constructor.
   655     /// The constructor of the class.
   590     ///
   656     ///
   591     /// \param graph The digraph the algorithm runs on.
   657     /// \param graph The digraph the algorithm runs on.
   592     NetworkSimplex(const GR& graph) :
   658     NetworkSimplex(const GR& graph) :
   593       _graph(graph),
   659       _graph(graph),
   594       _plower(NULL), _pupper(NULL), _pcost(NULL),
   660       _plower(NULL), _pupper(NULL), _pcost(NULL),
   595       _psupply(NULL), _pstsup(false),
   661       _psupply(NULL), _pstsup(false), _ptype(GEQ),
   596       _flow_map(NULL), _potential_map(NULL),
   662       _flow_map(NULL), _potential_map(NULL),
   597       _local_flow(false), _local_potential(false),
   663       _local_flow(false), _local_potential(false),
   598       _node_id(graph)
   664       _node_id(graph)
   599     {
   665     {
   600       LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
   666       LEMON_ASSERT(std::numeric_limits<Flow>::is_integer &&
   609     ~NetworkSimplex() {
   675     ~NetworkSimplex() {
   610       if (_local_flow) delete _flow_map;
   676       if (_local_flow) delete _flow_map;
   611       if (_local_potential) delete _potential_map;
   677       if (_local_potential) delete _potential_map;
   612     }
   678     }
   613 
   679 
       
   680     /// \name Parameters
       
   681     /// The parameters of the algorithm can be specified using these
       
   682     /// functions.
       
   683 
       
   684     /// @{
       
   685 
   614     /// \brief Set the lower bounds on the arcs.
   686     /// \brief Set the lower bounds on the arcs.
   615     ///
   687     ///
   616     /// This function sets the lower bounds on the arcs.
   688     /// This function sets the lower bounds on the arcs.
   617     /// If neither this function nor \ref boundMaps() is used before
   689     /// If neither this function nor \ref boundMaps() is used before
   618     /// calling \ref run(), the lower bounds will be set to zero
   690     /// calling \ref run(), the lower bounds will be set to zero
   758       _psource = s;
   830       _psource = s;
   759       _ptarget = t;
   831       _ptarget = t;
   760       _pstflow = k;
   832       _pstflow = k;
   761       return *this;
   833       return *this;
   762     }
   834     }
       
   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     }
   763 
   849 
   764     /// \brief Set the flow map.
   850     /// \brief Set the flow map.
   765     ///
   851     ///
   766     /// This function sets the flow map.
   852     /// This function sets the flow map.
   767     /// If it is not used before calling \ref run(), an instance will
   853     /// If it is not used before calling \ref run(), an instance will
   793         _local_potential = false;
   879         _local_potential = false;
   794       }
   880       }
   795       _potential_map = &map;
   881       _potential_map = &map;
   796       return *this;
   882       return *this;
   797     }
   883     }
       
   884     
       
   885     /// @}
   798 
   886 
   799     /// \name Execution Control
   887     /// \name Execution Control
   800     /// The algorithm can be executed using \ref run().
   888     /// The algorithm can be executed using \ref run().
   801 
   889 
   802     /// @{
   890     /// @{
   803 
   891 
   804     /// \brief Run the algorithm.
   892     /// \brief Run the algorithm.
   805     ///
   893     ///
   806     /// This function runs the algorithm.
   894     /// This function runs the algorithm.
   807     /// The paramters can be specified using \ref lowerMap(),
   895     /// The paramters can be specified using functions \ref lowerMap(),
   808     /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
   896     /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(),
   809     /// \ref costMap(), \ref supplyMap() and \ref stSupply()
   897     /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 
   810     /// functions. For example,
   898     /// \ref problemType(), \ref flowMap() and \ref potentialMap().
       
   899     /// For example,
   811     /// \code
   900     /// \code
   812     ///   NetworkSimplex<ListDigraph> ns(graph);
   901     ///   NetworkSimplex<ListDigraph> ns(graph);
   813     ///   ns.boundMaps(lower, upper).costMap(cost)
   902     ///   ns.boundMaps(lower, upper).costMap(cost)
   814     ///     .supplyMap(sup).run();
   903     ///     .supplyMap(sup).run();
   815     /// \endcode
   904     /// \endcode
   828     }
   917     }
   829 
   918 
   830     /// \brief Reset all the parameters that have been given before.
   919     /// \brief Reset all the parameters that have been given before.
   831     ///
   920     ///
   832     /// This function resets all the paramaters that have been given
   921     /// This function resets all the paramaters that have been given
   833     /// using \ref lowerMap(), \ref upperMap(), \ref capacityMap(),
   922     /// before using functions \ref lowerMap(), \ref upperMap(),
   834     /// \ref boundMaps(), \ref costMap(), \ref supplyMap() and
   923     /// \ref capacityMap(), \ref boundMaps(), \ref costMap(),
   835     /// \ref stSupply() functions before.
   924     /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 
       
   925     /// \ref flowMap() and \ref potentialMap().
   836     ///
   926     ///
   837     /// It is useful for multiple run() calls. If this function is not
   927     /// It is useful for multiple run() calls. If this function is not
   838     /// used, all the parameters given before are kept for the next
   928     /// used, all the parameters given before are kept for the next
   839     /// \ref run() call.
   929     /// \ref run() call.
   840     ///
   930     ///
   867       _plower = NULL;
   957       _plower = NULL;
   868       _pupper = NULL;
   958       _pupper = NULL;
   869       _pcost = NULL;
   959       _pcost = NULL;
   870       _psupply = NULL;
   960       _psupply = NULL;
   871       _pstsup = false;
   961       _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 
   872       return *this;
   970       return *this;
   873     }
   971     }
   874 
   972 
   875     /// @}
   973     /// @}
   876 
   974 
   998       _last_succ.resize(all_node_num);
  1096       _last_succ.resize(all_node_num);
   999       _state.resize(all_arc_num);
  1097       _state.resize(all_arc_num);
  1000 
  1098 
  1001       // Initialize node related data
  1099       // Initialize node related data
  1002       bool valid_supply = true;
  1100       bool valid_supply = true;
       
  1101       Flow sum_supply = 0;
  1003       if (!_pstsup && !_psupply) {
  1102       if (!_pstsup && !_psupply) {
  1004         _pstsup = true;
  1103         _pstsup = true;
  1005         _psource = _ptarget = NodeIt(_graph);
  1104         _psource = _ptarget = NodeIt(_graph);
  1006         _pstflow = 0;
  1105         _pstflow = 0;
  1007       }
  1106       }
  1008       if (_psupply) {
  1107       if (_psupply) {
  1009         Flow sum = 0;
       
  1010         int i = 0;
  1108         int i = 0;
  1011         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
  1109         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
  1012           _node_id[n] = i;
  1110           _node_id[n] = i;
  1013           _supply[i] = (*_psupply)[n];
  1111           _supply[i] = (*_psupply)[n];
  1014           sum += _supply[i];
  1112           sum_supply += _supply[i];
  1015         }
  1113         }
  1016         valid_supply = (sum == 0);
  1114         valid_supply = (_ptype == GEQ && sum_supply <= 0) ||
       
  1115                        (_ptype == LEQ && sum_supply >= 0);
  1017       } else {
  1116       } else {
  1018         int i = 0;
  1117         int i = 0;
  1019         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
  1118         for (NodeIt n(_graph); n != INVALID; ++n, ++i) {
  1020           _node_id[n] = i;
  1119           _node_id[n] = i;
  1021           _supply[i] = 0;
  1120           _supply[i] = 0;
  1022         }
  1121         }
  1023         _supply[_node_id[_psource]] =  _pstflow;
  1122         _supply[_node_id[_psource]] =  _pstflow;
  1024         _supply[_node_id[_ptarget]]   = -_pstflow;
  1123         _supply[_node_id[_ptarget]] = -_pstflow;
  1025       }
  1124       }
  1026       if (!valid_supply) return false;
  1125       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;
  1027 
  1211 
  1028       // Set data for the artificial root node
  1212       // Set data for the artificial root node
  1029       _root = _node_num;
  1213       _root = _node_num;
  1030       _parent[_root] = -1;
  1214       _parent[_root] = -1;
  1031       _pred[_root] = -1;
  1215       _pred[_root] = -1;
  1032       _thread[_root] = 0;
  1216       _thread[_root] = 0;
  1033       _rev_thread[0] = _root;
  1217       _rev_thread[0] = _root;
  1034       _succ_num[_root] = all_node_num;
  1218       _succ_num[_root] = all_node_num;
  1035       _last_succ[_root] = _root - 1;
  1219       _last_succ[_root] = _root - 1;
  1036       _supply[_root] = 0;
  1220       _supply[_root] = -sum_supply;
  1037       _pi[_root] = 0;
  1221       if (sum_supply < 0) {
       
  1222         _pi[_root] = -art_cost;
       
  1223       } else {
       
  1224         _pi[_root] = art_cost;
       
  1225       }
  1038 
  1226 
  1039       // Store the arcs in a mixed order
  1227       // Store the arcs in a mixed order
  1040       int k = std::max(int(sqrt(_arc_num)), 10);
  1228       int k = std::max(int(sqrt(_arc_num)), 10);
  1041       int i = 0;
  1229       int i = 0;
  1042       for (ArcIt e(_graph); e != INVALID; ++e) {
  1230       for (ArcIt e(_graph); e != INVALID; ++e) {
  1043         _arc_ref[i] = e;
  1231         _arc_ref[i] = e;
  1044         if ((i += k) >= _arc_num) i = (i % k) + 1;
  1232         if ((i += k) >= _arc_num) i = (i % k) + 1;
  1045       }
  1233       }
  1046 
  1234 
  1047       // Initialize arc maps
  1235       // 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();
       
  1052       if (_pupper && _pcost) {
  1236       if (_pupper && _pcost) {
  1053         for (int i = 0; i != _arc_num; ++i) {
  1237         for (int i = 0; i != _arc_num; ++i) {
  1054           Arc e = _arc_ref[i];
  1238           Arc e = _arc_ref[i];
  1055           _source[i] = _node_id[_graph.source(e)];
  1239           _source[i] = _node_id[_graph.source(e)];
  1056           _target[i] = _node_id[_graph.target(e)];
  1240           _target[i] = _node_id[_graph.target(e)];
  1081           for (int i = 0; i != _arc_num; ++i)
  1265           for (int i = 0; i != _arc_num; ++i)
  1082             _cost[i] = 1;
  1266             _cost[i] = 1;
  1083         }
  1267         }
  1084       }
  1268       }
  1085       
  1269       
  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 
       
  1098       // Remove non-zero lower bounds
  1270       // Remove non-zero lower bounds
  1099       if (_plower) {
  1271       if (_plower) {
  1100         for (int i = 0; i != _arc_num; ++i) {
  1272         for (int i = 0; i != _arc_num; ++i) {
  1101           Flow c = (*_plower)[_arc_ref[i]];
  1273           Flow c = (*_plower)[_arc_ref[i]];
  1102           if (c != 0) {
  1274           if (c != 0) {
  1116         _parent[u] = _root;
  1288         _parent[u] = _root;
  1117         _pred[u] = e;
  1289         _pred[u] = e;
  1118         _cost[e] = art_cost;
  1290         _cost[e] = art_cost;
  1119         _cap[e] = inf_cap;
  1291         _cap[e] = inf_cap;
  1120         _state[e] = STATE_TREE;
  1292         _state[e] = STATE_TREE;
  1121         if (_supply[u] >= 0) {
  1293         if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {
  1122           _flow[e] = _supply[u];
  1294           _flow[e] = _supply[u];
  1123           _forward[u] = true;
  1295           _forward[u] = true;
  1124           _pi[u] = -art_cost;
  1296           _pi[u] = -art_cost + _pi[_root];
  1125         } else {
  1297         } else {
  1126           _flow[e] = -_supply[u];
  1298           _flow[e] = -_supply[u];
  1127           _forward[u] = false;
  1299           _forward[u] = false;
  1128           _pi[u] = art_cost;
  1300           _pi[u] = art_cost + _pi[_root];
  1129         }
  1301         }
  1130       }
  1302       }
  1131 
  1303 
  1132       return true;
  1304       return true;
  1133     }
  1305     }
  1380           updateTreeStructure();
  1552           updateTreeStructure();
  1381           updatePotential();
  1553           updatePotential();
  1382         }
  1554         }
  1383       }
  1555       }
  1384 
  1556 
  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 
       
  1390       // Copy flow values to _flow_map
  1557       // Copy flow values to _flow_map
  1391       if (_plower) {
  1558       if (_plower) {
  1392         for (int i = 0; i != _arc_num; ++i) {
  1559         for (int i = 0; i != _arc_num; ++i) {
  1393           Arc e = _arc_ref[i];
  1560           Arc e = _arc_ref[i];
  1394           _flow_map->set(e, (*_plower)[e] + _flow[i]);
  1561           _flow_map->set(e, (*_plower)[e] + _flow[i]);