Changeset 687:6c408d864fa1 in lemon for lemon/network_simplex.h
 Timestamp:
 04/29/09 03:15:24 (11 years ago)
 Branch:
 default
 Phase:
 public
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/network_simplex.h
r665 r687 31 31 #include <lemon/core.h> 32 32 #include <lemon/math.h> 33 #include <lemon/maps.h>34 #include <lemon/circulation.h>35 #include <lemon/adaptors.h>36 33 37 34 namespace lemon { … … 51 48 /// In general this class is the fastest implementation available 52 49 /// 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 /// Moreover it supports both directions of the supply/demand inequality 51 /// constraints. For more information see \ref SupplyType. 52 /// 53 /// Most of the parameters of the problem (except for the digraph) 54 /// can be given using separate functions, and the algorithm can be 55 /// executed using the \ref run() function. If some parameters are not 56 /// specified, then default values will be used. 55 57 /// 56 58 /// \tparam GR The digraph type the algorithm runs on. … … 89 91 public: 90 92 91 /// \brief Enum type for selecting the pivot rule. 92 /// 93 /// Enum type for selecting the pivot rule for the \ref run() 94 /// function. 95 /// 96 /// \ref NetworkSimplex provides five different pivot rule 97 /// implementations that significantly affect the running time 98 /// of the algorithm. 99 /// By default \ref BLOCK_SEARCH "Block Search" is used, which 100 /// proved to be the most efficient and the most robust on various 101 /// test inputs according to our benchmark tests. 102 /// However another pivot rule can be selected using the \ref run() 103 /// function with the proper parameter. 104 enum PivotRule { 105 106 /// The First Eligible pivot rule. 107 /// The next eligible arc is selected in a wraparound fashion 108 /// in every iteration. 109 FIRST_ELIGIBLE, 110 111 /// The Best Eligible pivot rule. 112 /// The best eligible arc is selected in every iteration. 113 BEST_ELIGIBLE, 114 115 /// The Block Search pivot rule. 116 /// A specified number of arcs are examined in every iteration 117 /// in a wraparound fashion and the best eligible arc is selected 118 /// from this block. 119 BLOCK_SEARCH, 120 121 /// The Candidate List pivot rule. 122 /// In a major iteration a candidate list is built from eligible arcs 123 /// in a wraparound fashion and in the following minor iterations 124 /// the best eligible arc is selected from this list. 125 CANDIDATE_LIST, 126 127 /// The Altering Candidate List pivot rule. 128 /// It is a modified version of the Candidate List method. 129 /// It keeps only the several best eligible arcs from the former 130 /// candidate list and extends this list in every iteration. 131 ALTERING_LIST 93 /// \brief Problem type constants for the \c run() function. 94 /// 95 /// Enum type containing the problem type constants that can be 96 /// returned by the \ref run() function of the algorithm. 97 enum ProblemType { 98 /// The problem has no feasible solution (flow). 99 INFEASIBLE, 100 /// The problem has optimal solution (i.e. it is feasible and 101 /// bounded), and the algorithm has found optimal flow and node 102 /// potentials (primal and dual solutions). 103 OPTIMAL, 104 /// The objective function of the problem is unbounded, i.e. 105 /// there is a directed cycle having negative total cost and 106 /// infinite upper bound. 107 UNBOUNDED 132 108 }; 133 109 134 /// \brief Enum type for selecting the problem type.135 /// 136 /// Enum type for selecting the problem type, i.e. the direction of137 /// the inequalities in the supply/demand constraints of the138 /// \ref min_cost_flow "minimum cost flow problem".139 /// 140 /// The default problemtype is \c GEQ, since this form is supported110 /// \brief Constants for selecting the type of the supply constraints. 111 /// 112 /// Enum type containing constants for selecting the supply type, 113 /// i.e. the direction of the inequalities in the supply/demand 114 /// constraints of the \ref min_cost_flow "minimum cost flow problem". 115 /// 116 /// The default supply type is \c GEQ, since this form is supported 141 117 /// 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()118 /// algorithm, as well. 119 /// The \c LEQ problem type can be selected using the \ref supplyType() 144 120 /// function. 145 121 /// 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 the151 /// problem is the following.122 /// Note that the equality form is a special case of both supply types. 123 enum SupplyType { 124 125 /// This option means that there are <em>"greater or equal"</em> 126 /// supply/demand constraints in the definition, i.e. the exact 127 /// formulation of the problem is the following. 152 128 /** 153 129 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] … … 165 141 CARRY_SUPPLIES = GEQ, 166 142 167 /// This option means that there are "<em>less or equal</em>"168 /// constraints in the defintion, i.e. the exact formulation of the169 /// problem is the following.143 /// This option means that there are <em>"less or equal"</em> 144 /// supply/demand constraints in the definition, i.e. the exact 145 /// formulation of the problem is the following. 170 146 /** 171 147 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] … … 183 159 SATISFY_DEMANDS = LEQ 184 160 }; 185 161 162 /// \brief Constants for selecting the pivot rule. 163 /// 164 /// Enum type containing constants for selecting the pivot rule for 165 /// the \ref run() function. 166 /// 167 /// \ref NetworkSimplex provides five different pivot rule 168 /// implementations that significantly affect the running time 169 /// of the algorithm. 170 /// By default \ref BLOCK_SEARCH "Block Search" is used, which 171 /// proved to be the most efficient and the most robust on various 172 /// test inputs according to our benchmark tests. 173 /// However another pivot rule can be selected using the \ref run() 174 /// function with the proper parameter. 175 enum PivotRule { 176 177 /// The First Eligible pivot rule. 178 /// The next eligible arc is selected in a wraparound fashion 179 /// in every iteration. 180 FIRST_ELIGIBLE, 181 182 /// The Best Eligible pivot rule. 183 /// The best eligible arc is selected in every iteration. 184 BEST_ELIGIBLE, 185 186 /// The Block Search pivot rule. 187 /// A specified number of arcs are examined in every iteration 188 /// in a wraparound fashion and the best eligible arc is selected 189 /// from this block. 190 BLOCK_SEARCH, 191 192 /// The Candidate List pivot rule. 193 /// In a major iteration a candidate list is built from eligible arcs 194 /// in a wraparound fashion and in the following minor iterations 195 /// the best eligible arc is selected from this list. 196 CANDIDATE_LIST, 197 198 /// The Altering Candidate List pivot rule. 199 /// It is a modified version of the Candidate List method. 200 /// It keeps only the several best eligible arcs from the former 201 /// candidate list and extends this list in every iteration. 202 ALTERING_LIST 203 }; 204 186 205 private: 187 206 … … 221 240 Node _psource, _ptarget; 222 241 Flow _pstflow; 223 ProblemType _ptype; 242 SupplyType _stype; 243 244 Flow _sum_supply; 224 245 225 246 // Result maps … … 259 280 int stem, par_stem, new_stem; 260 281 Flow delta; 282 283 public: 284 285 /// \brief Constant for infinite upper bounds (capacities). 286 /// 287 /// Constant for infinite upper bounds (capacities). 288 /// It is \c std::numeric_limits<Flow>::infinity() if available, 289 /// \c std::numeric_limits<Flow>::max() otherwise. 290 const Flow INF; 261 291 262 292 private: … … 662 692 _graph(graph), 663 693 _plower(NULL), _pupper(NULL), _pcost(NULL), 664 _psupply(NULL), _pstsup(false), _ ptype(GEQ),694 _psupply(NULL), _pstsup(false), _stype(GEQ), 665 695 _flow_map(NULL), _potential_map(NULL), 666 696 _local_flow(false), _local_potential(false), 667 _node_id(graph) 697 _node_id(graph), 698 INF(std::numeric_limits<Flow>::has_infinity ? 699 std::numeric_limits<Flow>::infinity() : 700 std::numeric_limits<Flow>::max()) 668 701 { 669 LEMON_ASSERT(std::numeric_limits<Flow>::is_integer && 670 std::numeric_limits<Flow>::is_signed, 671 "The flow type of NetworkSimplex must be signed integer"); 672 LEMON_ASSERT(std::numeric_limits<Cost>::is_integer && 673 std::numeric_limits<Cost>::is_signed, 674 "The cost type of NetworkSimplex must be signed integer"); 702 // Check the value types 703 LEMON_ASSERT(std::numeric_limits<Flow>::is_signed, 704 "The flow type of NetworkSimplex must be signed"); 705 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, 706 "The cost type of NetworkSimplex must be signed"); 675 707 } 676 708 … … 690 722 /// 691 723 /// This function sets the lower bounds on the arcs. 692 /// If neither this function nor \ref boundMaps() is used before 693 /// calling \ref run(), the lower bounds will be set to zero 694 /// on all arcs. 724 /// If it is not used before calling \ref run(), the lower bounds 725 /// will be set to zero on all arcs. 695 726 /// 696 727 /// \param map An arc map storing the lower bounds. … … 699 730 /// 700 731 /// \return <tt>(*this)</tt> 701 template <typename L OWER>702 NetworkSimplex& lowerMap(const L OWER& map) {732 template <typename LowerMap> 733 NetworkSimplex& lowerMap(const LowerMap& map) { 703 734 delete _plower; 704 735 _plower = new FlowArcMap(_graph); … … 712 743 /// 713 744 /// This function sets the upper bounds (capacities) on the arcs. 714 /// If none of the functions \ref upperMap(), \ref capacityMap() 715 /// and \ref boundMaps() is used before calling \ref run(), 716 /// the upper bounds (capacities) will be set to 717 /// \c std::numeric_limits<Flow>::max() on all arcs. 745 /// If it is not used before calling \ref run(), the upper bounds 746 /// will be set to \ref INF on all arcs (i.e. the flow value will be 747 /// unbounded from above on each arc). 718 748 /// 719 749 /// \param map An arc map storing the upper bounds. … … 722 752 /// 723 753 /// \return <tt>(*this)</tt> 724 template<typename U PPER>725 NetworkSimplex& upperMap(const U PPER& map) {754 template<typename UpperMap> 755 NetworkSimplex& upperMap(const UpperMap& map) { 726 756 delete _pupper; 727 757 _pupper = new FlowArcMap(_graph); … … 732 762 } 733 763 734 /// \brief Set the upper bounds (capacities) on the arcs.735 ///736 /// This function sets the upper bounds (capacities) on the arcs.737 /// It is just an alias for \ref upperMap().738 ///739 /// \return <tt>(*this)</tt>740 template<typename CAP>741 NetworkSimplex& capacityMap(const CAP& map) {742 return upperMap(map);743 }744 745 /// \brief Set the lower and upper bounds on the arcs.746 ///747 /// This function sets the lower and upper bounds on the arcs.748 /// If neither this function nor \ref lowerMap() is used before749 /// calling \ref run(), the lower bounds will be set to zero750 /// on all arcs.751 /// If none of the functions \ref upperMap(), \ref capacityMap()752 /// and \ref boundMaps() is used before calling \ref run(),753 /// the upper bounds (capacities) will be set to754 /// \c std::numeric_limits<Flow>::max() on all arcs.755 ///756 /// \param lower An arc map storing the lower bounds.757 /// \param upper An arc map storing the upper bounds.758 ///759 /// The \c Value type of the maps must be convertible to the760 /// \c Flow type of the algorithm.761 ///762 /// \note This function is just a shortcut of calling \ref lowerMap()763 /// and \ref upperMap() separately.764 ///765 /// \return <tt>(*this)</tt>766 template <typename LOWER, typename UPPER>767 NetworkSimplex& boundMaps(const LOWER& lower, const UPPER& upper) {768 return lowerMap(lower).upperMap(upper);769 }770 771 764 /// \brief Set the costs of the arcs. 772 765 /// … … 780 773 /// 781 774 /// \return <tt>(*this)</tt> 782 template<typename C OST>783 NetworkSimplex& costMap(const C OST& map) {775 template<typename CostMap> 776 NetworkSimplex& costMap(const CostMap& map) { 784 777 delete _pcost; 785 778 _pcost = new CostArcMap(_graph); … … 802 795 /// 803 796 /// \return <tt>(*this)</tt> 804 template<typename S UP>805 NetworkSimplex& supplyMap(const S UP& map) {797 template<typename SupplyMap> 798 NetworkSimplex& supplyMap(const SupplyMap& map) { 806 799 delete _psupply; 807 800 _pstsup = false; … … 820 813 /// calling \ref run(), the supply of each node will be set to zero. 821 814 /// (It makes sense only if nonzero lower bounds are given.) 815 /// 816 /// Using this function has the same effect as using \ref supplyMap() 817 /// with such a map in which \c k is assigned to \c s, \c k is 818 /// assigned to \c t and all other nodes have zero supply value. 822 819 /// 823 820 /// \param s The source node. … … 837 834 } 838 835 839 /// \brief Set the problem type.840 /// 841 /// This function sets the problem type for the algorithm.842 /// If it is not used before calling \ref run(), the \ref GEQ problem836 /// \brief Set the type of the supply constraints. 837 /// 838 /// This function sets the type of the supply/demand constraints. 839 /// If it is not used before calling \ref run(), the \ref GEQ supply 843 840 /// type will be used. 844 841 /// 845 /// For more information see \ref ProblemType.842 /// For more information see \ref SupplyType. 846 843 /// 847 844 /// \return <tt>(*this)</tt> 848 NetworkSimplex& problemType(ProblemType problem_type) {849 _ ptype = problem_type;845 NetworkSimplex& supplyType(SupplyType supply_type) { 846 _stype = supply_type; 850 847 return *this; 851 848 } … … 897 894 /// This function runs the algorithm. 898 895 /// The paramters can be specified using functions \ref lowerMap(), 899 /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(), 900 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 901 /// \ref problemType(), \ref flowMap() and \ref potentialMap(). 896 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 897 /// \ref supplyType(), \ref flowMap() and \ref potentialMap(). 902 898 /// For example, 903 899 /// \code 904 900 /// NetworkSimplex<ListDigraph> ns(graph); 905 /// ns. boundMaps(lower,upper).costMap(cost)901 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 906 902 /// .supplyMap(sup).run(); 907 903 /// \endcode … … 915 911 /// algorithm. For more information see \ref PivotRule. 916 912 /// 917 /// \return \c true if a feasible flow can be found. 918 bool run(PivotRule pivot_rule = BLOCK_SEARCH) { 919 return init() && start(pivot_rule); 913 /// \return \c INFEASIBLE if no feasible flow exists, 914 /// \n \c OPTIMAL if the problem has optimal solution 915 /// (i.e. it is feasible and bounded), and the algorithm has found 916 /// optimal flow and node potentials (primal and dual solutions), 917 /// \n \c UNBOUNDED if the objective function of the problem is 918 /// unbounded, i.e. there is a directed cycle having negative total 919 /// cost and infinite upper bound. 920 /// 921 /// \see ProblemType, PivotRule 922 ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { 923 if (!init()) return INFEASIBLE; 924 return start(pivot_rule); 920 925 } 921 926 … … 924 929 /// This function resets all the paramaters that have been given 925 930 /// before using functions \ref lowerMap(), \ref upperMap(), 926 /// \ref capacityMap(), \ref boundMaps(), \ref costMap(), 927 /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 931 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(), 928 932 /// \ref flowMap() and \ref potentialMap(). 929 933 /// … … 937 941 /// 938 942 /// // First run 939 /// ns.lowerMap(lower). capacityMap(cap).costMap(cost)943 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 940 944 /// .supplyMap(sup).run(); 941 945 /// … … 948 952 /// // (the lower bounds will be set to zero on all arcs) 949 953 /// ns.reset(); 950 /// ns. capacityMap(cap).costMap(cost)954 /// ns.upperMap(capacity).costMap(cost) 951 955 /// .supplyMap(sup).run(); 952 956 /// \endcode … … 963 967 _psupply = NULL; 964 968 _pstsup = false; 965 _ ptype = GEQ;969 _stype = GEQ; 966 970 if (_local_flow) delete _flow_map; 967 971 if (_local_potential) delete _potential_map; … … 986 990 /// 987 991 /// This function returns the total cost of the found flow. 988 /// The complexity of the functionis O(e).992 /// Its complexity is O(e). 989 993 /// 990 994 /// \note The return type of the function can be specified as a … … 998 1002 /// 999 1003 /// \pre \ref run() must be called before using this function. 1000 template <typename Num>1001 NumtotalCost() const {1002 Numc = 0;1004 template <typename Value> 1005 Value totalCost() const { 1006 Value c = 0; 1003 1007 if (_pcost) { 1004 1008 for (ArcIt e(_graph); e != INVALID; ++e) … … 1051 1055 /// This function returns a const reference to a node map storing 1052 1056 /// the found potentials, which form the dual solution of the 1053 /// \ref min_cost_flow "minimum cost flow " problem.1057 /// \ref min_cost_flow "minimum cost flow problem". 1054 1058 /// 1055 1059 /// \pre \ref run() must be called before using this function. … … 1102 1106 // Initialize node related data 1103 1107 bool valid_supply = true; 1104 Flowsum_supply = 0;1108 _sum_supply = 0; 1105 1109 if (!_pstsup && !_psupply) { 1106 1110 _pstsup = true; … … 1113 1117 _node_id[n] = i; 1114 1118 _supply[i] = (*_psupply)[n]; 1115 sum_supply += _supply[i];1116 } 1117 valid_supply = (_ ptype == GEQ &&sum_supply <= 0) 1118 (_ ptype == LEQ &&sum_supply >= 0);1119 _sum_supply += _supply[i]; 1120 } 1121 valid_supply = (_stype == GEQ && _sum_supply <= 0)  1122 (_stype == LEQ && _sum_supply >= 0); 1119 1123 } else { 1120 1124 int i = 0; … … 1128 1132 if (!valid_supply) return false; 1129 1133 1130 // Infinite capacity value1131 Flow inf_cap =1132 std::numeric_limits<Flow>::has_infinity ?1133 std::numeric_limits<Flow>::infinity() :1134 std::numeric_limits<Flow>::max();1135 1136 1134 // Initialize artifical cost 1137 Cost art_cost;1135 Cost ART_COST; 1138 1136 if (std::numeric_limits<Cost>::is_exact) { 1139 art_cost= std::numeric_limits<Cost>::max() / 4 + 1;1137 ART_COST = std::numeric_limits<Cost>::max() / 4 + 1; 1140 1138 } else { 1141 art_cost= std::numeric_limits<Cost>::min();1139 ART_COST = std::numeric_limits<Cost>::min(); 1142 1140 for (int i = 0; i != _arc_num; ++i) { 1143 if (_cost[i] > art_cost) art_cost = _cost[i]; 1144 } 1145 art_cost = (art_cost + 1) * _node_num; 1146 } 1147 1148 // Run Circulation to check if a feasible solution exists 1149 typedef ConstMap<Arc, Flow> ConstArcMap; 1150 ConstArcMap zero_arc_map(0), inf_arc_map(inf_cap); 1151 FlowNodeMap *csup = NULL; 1152 bool local_csup = false; 1153 if (_psupply) { 1154 csup = _psupply; 1155 } else { 1156 csup = new FlowNodeMap(_graph, 0); 1157 (*csup)[_psource] = _pstflow; 1158 (*csup)[_ptarget] = _pstflow; 1159 local_csup = true; 1160 } 1161 bool circ_result = false; 1162 if (_ptype == GEQ  (_ptype == LEQ && sum_supply == 0)) { 1163 // GEQ problem type 1164 if (_plower) { 1165 if (_pupper) { 1166 Circulation<GR, FlowArcMap, FlowArcMap, FlowNodeMap> 1167 circ(_graph, *_plower, *_pupper, *csup); 1168 circ_result = circ.run(); 1169 } else { 1170 Circulation<GR, FlowArcMap, ConstArcMap, FlowNodeMap> 1171 circ(_graph, *_plower, inf_arc_map, *csup); 1172 circ_result = circ.run(); 1173 } 1174 } else { 1175 if (_pupper) { 1176 Circulation<GR, ConstArcMap, FlowArcMap, FlowNodeMap> 1177 circ(_graph, zero_arc_map, *_pupper, *csup); 1178 circ_result = circ.run(); 1179 } else { 1180 Circulation<GR, ConstArcMap, ConstArcMap, FlowNodeMap> 1181 circ(_graph, zero_arc_map, inf_arc_map, *csup); 1182 circ_result = circ.run(); 1183 } 1184 } 1185 } else { 1186 // LEQ problem type 1187 typedef ReverseDigraph<const GR> RevGraph; 1188 typedef NegMap<FlowNodeMap> NegNodeMap; 1189 RevGraph rgraph(_graph); 1190 NegNodeMap neg_csup(*csup); 1191 if (_plower) { 1192 if (_pupper) { 1193 Circulation<RevGraph, FlowArcMap, FlowArcMap, NegNodeMap> 1194 circ(rgraph, *_plower, *_pupper, neg_csup); 1195 circ_result = circ.run(); 1196 } else { 1197 Circulation<RevGraph, FlowArcMap, ConstArcMap, NegNodeMap> 1198 circ(rgraph, *_plower, inf_arc_map, neg_csup); 1199 circ_result = circ.run(); 1200 } 1201 } else { 1202 if (_pupper) { 1203 Circulation<RevGraph, ConstArcMap, FlowArcMap, NegNodeMap> 1204 circ(rgraph, zero_arc_map, *_pupper, neg_csup); 1205 circ_result = circ.run(); 1206 } else { 1207 Circulation<RevGraph, ConstArcMap, ConstArcMap, NegNodeMap> 1208 circ(rgraph, zero_arc_map, inf_arc_map, neg_csup); 1209 circ_result = circ.run(); 1210 } 1211 } 1212 } 1213 if (local_csup) delete csup; 1214 if (!circ_result) return false; 1141 if (_cost[i] > ART_COST) ART_COST = _cost[i]; 1142 } 1143 ART_COST = (ART_COST + 1) * _node_num; 1144 } 1215 1145 1216 1146 // Set data for the artificial root node … … 1222 1152 _succ_num[_root] = all_node_num; 1223 1153 _last_succ[_root] = _root  1; 1224 _supply[_root] =  sum_supply;1225 if ( sum_supply < 0) {1226 _pi[_root] =  art_cost;1154 _supply[_root] = _sum_supply; 1155 if (_sum_supply < 0) { 1156 _pi[_root] = ART_COST; 1227 1157 } else { 1228 _pi[_root] = art_cost;1158 _pi[_root] = ART_COST; 1229 1159 } 1230 1160 … … 1261 1191 } else { 1262 1192 for (int i = 0; i != _arc_num; ++i) 1263 _cap[i] = inf_cap;1193 _cap[i] = INF; 1264 1194 } 1265 1195 if (_pcost) { … … 1276 1206 for (int i = 0; i != _arc_num; ++i) { 1277 1207 Flow c = (*_plower)[_arc_ref[i]]; 1278 if (c != 0) { 1279 _cap[i] = c; 1208 if (c > 0) { 1209 if (_cap[i] < INF) _cap[i] = c; 1210 _supply[_source[i]] = c; 1211 _supply[_target[i]] += c; 1212 } 1213 else if (c < 0) { 1214 if (_cap[i] < INF + c) { 1215 _cap[i] = c; 1216 } else { 1217 _cap[i] = INF; 1218 } 1280 1219 _supply[_source[i]] = c; 1281 1220 _supply[_target[i]] += c; … … 1292 1231 _parent[u] = _root; 1293 1232 _pred[u] = e; 1294 _cost[e] = art_cost;1295 _cap[e] = inf_cap;1233 _cost[e] = ART_COST; 1234 _cap[e] = INF; 1296 1235 _state[e] = STATE_TREE; 1297 if (_supply[u] > 0  (_supply[u] == 0 && sum_supply <= 0)) {1236 if (_supply[u] > 0  (_supply[u] == 0 && _sum_supply <= 0)) { 1298 1237 _flow[e] = _supply[u]; 1299 1238 _forward[u] = true; 1300 _pi[u] =  art_cost+ _pi[_root];1239 _pi[u] = ART_COST + _pi[_root]; 1301 1240 } else { 1302 1241 _flow[e] = _supply[u]; 1303 1242 _forward[u] = false; 1304 _pi[u] = art_cost+ _pi[_root];1243 _pi[u] = ART_COST + _pi[_root]; 1305 1244 } 1306 1245 } … … 1343 1282 for (int u = first; u != join; u = _parent[u]) { 1344 1283 e = _pred[u]; 1345 d = _forward[u] ? _flow[e] : _cap[e]  _flow[e]; 1284 d = _forward[u] ? 1285 _flow[e] : (_cap[e] == INF ? INF : _cap[e]  _flow[e]); 1346 1286 if (d < delta) { 1347 1287 delta = d; … … 1353 1293 for (int u = second; u != join; u = _parent[u]) { 1354 1294 e = _pred[u]; 1355 d = _forward[u] ? _cap[e]  _flow[e] : _flow[e]; 1295 d = _forward[u] ? 1296 (_cap[e] == INF ? INF : _cap[e]  _flow[e]) : _flow[e]; 1356 1297 if (d <= delta) { 1357 1298 delta = d; … … 1527 1468 1528 1469 // Execute the algorithm 1529 boolstart(PivotRule pivot_rule) {1470 ProblemType start(PivotRule pivot_rule) { 1530 1471 // Select the pivot rule implementation 1531 1472 switch (pivot_rule) { … … 1541 1482 return start<AlteringListPivotRule>(); 1542 1483 } 1543 return false;1484 return INFEASIBLE; // avoid warning 1544 1485 } 1545 1486 1546 1487 template <typename PivotRuleImpl> 1547 boolstart() {1488 ProblemType start() { 1548 1489 PivotRuleImpl pivot(*this); 1549 1490 … … 1552 1493 findJoinNode(); 1553 1494 bool change = findLeavingArc(); 1495 if (delta >= INF) return UNBOUNDED; 1554 1496 changeFlow(change); 1555 1497 if (change) { 1556 1498 updateTreeStructure(); 1557 1499 updatePotential(); 1500 } 1501 } 1502 1503 // Check feasibility 1504 if (_sum_supply < 0) { 1505 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1506 if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE; 1507 } 1508 } 1509 else if (_sum_supply > 0) { 1510 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1511 if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE; 1512 } 1513 } 1514 else { 1515 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1516 if (_flow[e] != 0) return INFEASIBLE; 1558 1517 } 1559 1518 } … … 1575 1534 } 1576 1535 1577 return true;1536 return OPTIMAL; 1578 1537 } 1579 1538
Note: See TracChangeset
for help on using the changeset viewer.