Changeset 656:e6927fe719e6 in lemon for lemon
- Timestamp:
- 04/17/09 18:04:36 (16 years ago)
- Branch:
- default
- Phase:
- public
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/network_simplex.h
r655 r656 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> 33 36 34 37 namespace lemon { … … 48 51 /// In general this class is the fastest implementation available 49 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 56 /// \tparam GR The digraph type the algorithm runs on. … … 59 64 /// 60 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 68 template <typename GR, typename F = int, typename C = F> 63 69 class NetworkSimplex … … 69 75 /// The cost type of the algorithm 70 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 83 /// The type of the flow map 72 84 typedef typename GR::template ArcMap<Flow> FlowMap; 73 85 /// The type of the potential map 74 86 typedef typename GR::template NodeMap<Cost> PotentialMap; 87 #endif 75 88 76 89 public: … … 118 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 186 private: … … 156 221 Node _psource, _ptarget; 157 222 Flow _pstflow; 223 ProblemType _ptype; 158 224 159 225 // Result maps … … 587 653 /// \brief Constructor. 588 654 /// 589 /// Constructor.655 /// The constructor of the class. 590 656 /// 591 657 /// \param graph The digraph the algorithm runs on. … … 593 659 _graph(graph), 594 660 _plower(NULL), _pupper(NULL), _pcost(NULL), 595 _psupply(NULL), _pstsup(false), 661 _psupply(NULL), _pstsup(false), _ptype(GEQ), 596 662 _flow_map(NULL), _potential_map(NULL), 597 663 _local_flow(false), _local_potential(false), … … 612 678 } 613 679 680 /// \name Parameters 681 /// The parameters of the algorithm can be specified using these 682 /// functions. 683 684 /// @{ 685 614 686 /// \brief Set the lower bounds on the arcs. 615 687 /// … … 761 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 850 /// \brief Set the flow map. … … 796 882 return *this; 797 883 } 884 885 /// @} 798 886 799 887 /// \name Execution Control … … 805 893 /// 806 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 896 /// \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, 811 900 /// \code 812 901 /// NetworkSimplex<ListDigraph> ns(graph); … … 831 920 /// 832 921 /// 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(). 836 926 /// 837 927 /// It is useful for multiple run() calls. If this function is not … … 870 960 _psupply = NULL; 871 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 970 return *this; 873 971 } … … 1001 1099 // Initialize node related data 1002 1100 bool valid_supply = true; 1101 Flow sum_supply = 0; 1003 1102 if (!_pstsup && !_psupply) { 1004 1103 _pstsup = true; … … 1007 1106 } 1008 1107 if (_psupply) { 1009 Flow sum = 0;1010 1108 int i = 0; 1011 1109 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 1012 1110 _node_id[n] = i; 1013 1111 _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); 1017 1116 } else { 1018 1117 int i = 0; … … 1022 1121 } 1023 1122 _supply[_node_id[_psource]] = _pstflow; 1024 _supply[_node_id[_ptarget]] 1123 _supply[_node_id[_ptarget]] = -_pstflow; 1025 1124 } 1026 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 1212 // Set data for the artificial root node … … 1034 1218 _succ_num[_root] = all_node_num; 1035 1219 _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 } 1038 1226 1039 1227 // Store the arcs in a mixed order … … 1046 1234 1047 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 1236 if (_pupper && _pcost) { 1053 1237 for (int i = 0; i != _arc_num; ++i) { … … 1084 1268 } 1085 1269 1086 // Initialize artifical cost1087 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 1270 // Remove non-zero lower bounds 1099 1271 if (_plower) { … … 1119 1291 _cap[e] = inf_cap; 1120 1292 _state[e] = STATE_TREE; 1121 if (_supply[u] > = 0) {1293 if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) { 1122 1294 _flow[e] = _supply[u]; 1123 1295 _forward[u] = true; 1124 _pi[u] = -art_cost ;1296 _pi[u] = -art_cost + _pi[_root]; 1125 1297 } else { 1126 1298 _flow[e] = -_supply[u]; 1127 1299 _forward[u] = false; 1128 _pi[u] = art_cost ;1300 _pi[u] = art_cost + _pi[_root]; 1129 1301 } 1130 1302 } … … 1383 1555 } 1384 1556 1385 // Check if the flow amount equals zero on all the artificial arcs1386 for (int e = _arc_num; e != _arc_num + _node_num; ++e) {1387 if (_flow[e] > 0) return false;1388 }1389 1390 1557 // Copy flow values to _flow_map 1391 1558 if (_plower) {
Note: See TracChangeset
for help on using the changeset viewer.