Changeset 689:111698359429 in lemon for lemon/network_simplex.h
 Timestamp:
 04/29/09 16:54:27 (10 years ago)
 Branch:
 default
 Phase:
 public
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

lemon/network_simplex.h
r688 r689 73 73 public: 74 74 75 /// The flow type of the algorithm75 /// The type of the flow amounts, capacity bounds and supply values 76 76 typedef V Value; 77 /// The cost type of the algorithm77 /// The type of the arc costs 78 78 typedef C Cost; 79 #ifdef DOXYGEN80 /// The type of the flow map81 typedef GR::ArcMap<Value> FlowMap;82 /// The type of the potential map83 typedef GR::NodeMap<Cost> PotentialMap;84 #else85 /// The type of the flow map86 typedef typename GR::template ArcMap<Value> FlowMap;87 /// The type of the potential map88 typedef typename GR::template NodeMap<Cost> PotentialMap;89 #endif90 79 91 80 public: … … 207 196 TEMPLATE_DIGRAPH_TYPEDEFS(GR); 208 197 209 typedef typename GR::template ArcMap<Value> ValueArcMap;210 typedef typename GR::template ArcMap<Cost> CostArcMap;211 typedef typename GR::template NodeMap<Value> ValueNodeMap;212 213 198 typedef std::vector<Arc> ArcVector; 214 199 typedef std::vector<Node> NodeVector; 215 200 typedef std::vector<int> IntVector; 216 201 typedef std::vector<bool> BoolVector; 217 typedef std::vector<Value> FlowVector;202 typedef std::vector<Value> ValueVector; 218 203 typedef std::vector<Cost> CostVector; 219 204 … … 233 218 234 219 // Parameters of the problem 235 ValueArcMap *_plower; 236 ValueArcMap *_pupper; 237 CostArcMap *_pcost; 238 ValueNodeMap *_psupply; 239 bool _pstsup; 240 Node _psource, _ptarget; 241 Value _pstflow; 220 bool _have_lower; 242 221 SupplyType _stype; 243 244 222 Value _sum_supply; 245 246 // Result maps247 FlowMap *_flow_map;248 PotentialMap *_potential_map;249 bool _local_flow;250 bool _local_potential;251 223 252 224 // Data structures for storing the digraph 253 225 IntNodeMap _node_id; 254 ArcVector _arc_ref;226 IntArcMap _arc_id; 255 227 IntVector _source; 256 228 IntVector _target; 257 229 258 230 // Node and arc data 259 FlowVector _cap; 231 ValueVector _lower; 232 ValueVector _upper; 233 ValueVector _cap; 260 234 CostVector _cost; 261 FlowVector _supply;262 FlowVector _flow;235 ValueVector _supply; 236 ValueVector _flow; 263 237 CostVector _pi; 264 238 … … 690 664 /// \param graph The digraph the algorithm runs on. 691 665 NetworkSimplex(const GR& graph) : 692 _graph(graph), 693 _plower(NULL), _pupper(NULL), _pcost(NULL), 694 _psupply(NULL), _pstsup(false), _stype(GEQ), 695 _flow_map(NULL), _potential_map(NULL), 696 _local_flow(false), _local_potential(false), 697 _node_id(graph), 666 _graph(graph), _node_id(graph), _arc_id(graph), 698 667 INF(std::numeric_limits<Value>::has_infinity ? 699 668 std::numeric_limits<Value>::infinity() : … … 705 674 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, 706 675 "The cost type of NetworkSimplex must be signed"); 707 } 708 709 /// Destructor. 710 ~NetworkSimplex() { 711 if (_local_flow) delete _flow_map; 712 if (_local_potential) delete _potential_map; 676 677 // Resize vectors 678 _node_num = countNodes(_graph); 679 _arc_num = countArcs(_graph); 680 int all_node_num = _node_num + 1; 681 int all_arc_num = _arc_num + _node_num; 682 683 _source.resize(all_arc_num); 684 _target.resize(all_arc_num); 685 686 _lower.resize(all_arc_num); 687 _upper.resize(all_arc_num); 688 _cap.resize(all_arc_num); 689 _cost.resize(all_arc_num); 690 _supply.resize(all_node_num); 691 _flow.resize(all_arc_num); 692 _pi.resize(all_node_num); 693 694 _parent.resize(all_node_num); 695 _pred.resize(all_node_num); 696 _forward.resize(all_node_num); 697 _thread.resize(all_node_num); 698 _rev_thread.resize(all_node_num); 699 _succ_num.resize(all_node_num); 700 _last_succ.resize(all_node_num); 701 _state.resize(all_arc_num); 702 703 // Copy the graph (store the arcs in a mixed order) 704 int i = 0; 705 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 706 _node_id[n] = i; 707 } 708 int k = std::max(int(std::sqrt(double(_arc_num))), 10); 709 i = 0; 710 for (ArcIt a(_graph); a != INVALID; ++a) { 711 _arc_id[a] = i; 712 _source[i] = _node_id[_graph.source(a)]; 713 _target[i] = _node_id[_graph.target(a)]; 714 if ((i += k) >= _arc_num) i = (i % k) + 1; 715 } 716 717 // Initialize maps 718 for (int i = 0; i != _node_num; ++i) { 719 _supply[i] = 0; 720 } 721 for (int i = 0; i != _arc_num; ++i) { 722 _lower[i] = 0; 723 _upper[i] = INF; 724 _cost[i] = 1; 725 } 726 _have_lower = false; 727 _stype = GEQ; 713 728 } 714 729 … … 732 747 template <typename LowerMap> 733 748 NetworkSimplex& lowerMap(const LowerMap& map) { 734 delete _plower; 735 _plower = new ValueArcMap(_graph); 749 _have_lower = true; 736 750 for (ArcIt a(_graph); a != INVALID; ++a) { 737 (*_plower)[a] = map[a];751 _lower[_arc_id[a]] = map[a]; 738 752 } 739 753 return *this; … … 754 768 template<typename UpperMap> 755 769 NetworkSimplex& upperMap(const UpperMap& map) { 756 delete _pupper;757 _pupper = new ValueArcMap(_graph);758 770 for (ArcIt a(_graph); a != INVALID; ++a) { 759 (*_pupper)[a] = map[a];771 _upper[_arc_id[a]] = map[a]; 760 772 } 761 773 return *this; … … 775 787 template<typename CostMap> 776 788 NetworkSimplex& costMap(const CostMap& map) { 777 delete _pcost;778 _pcost = new CostArcMap(_graph);779 789 for (ArcIt a(_graph); a != INVALID; ++a) { 780 (*_pcost)[a] = map[a];790 _cost[_arc_id[a]] = map[a]; 781 791 } 782 792 return *this; … … 797 807 template<typename SupplyMap> 798 808 NetworkSimplex& supplyMap(const SupplyMap& map) { 799 delete _psupply;800 _pstsup = false;801 _psupply = new ValueNodeMap(_graph);802 809 for (NodeIt n(_graph); n != INVALID; ++n) { 803 (*_psupply)[n] = map[n];810 _supply[_node_id[n]] = map[n]; 804 811 } 805 812 return *this; … … 825 832 /// \return <tt>(*this)</tt> 826 833 NetworkSimplex& stSupply(const Node& s, const Node& t, Value k) { 827 delete _psupply; 828 _psupply = NULL; 829 _pstsup = true; 830 _psource = s; 831 _ptarget = t; 832 _pstflow = k; 834 for (int i = 0; i != _node_num; ++i) { 835 _supply[i] = 0; 836 } 837 _supply[_node_id[s]] = k; 838 _supply[_node_id[t]] = k; 833 839 return *this; 834 840 } … … 848 854 } 849 855 850 /// \brief Set the flow map.851 ///852 /// This function sets the flow map.853 /// If it is not used before calling \ref run(), an instance will854 /// be allocated automatically. The destructor deallocates this855 /// automatically allocated map, of course.856 ///857 /// \return <tt>(*this)</tt>858 NetworkSimplex& flowMap(FlowMap& map) {859 if (_local_flow) {860 delete _flow_map;861 _local_flow = false;862 }863 _flow_map = ↦864 return *this;865 }866 867 /// \brief Set the potential map.868 ///869 /// This function sets the potential map, which is used for storing870 /// the dual solution.871 /// If it is not used before calling \ref run(), an instance will872 /// be allocated automatically. The destructor deallocates this873 /// automatically allocated map, of course.874 ///875 /// \return <tt>(*this)</tt>876 NetworkSimplex& potentialMap(PotentialMap& map) {877 if (_local_potential) {878 delete _potential_map;879 _local_potential = false;880 }881 _potential_map = ↦882 return *this;883 }884 885 856 /// @} 886 857 … … 895 866 /// The paramters can be specified using functions \ref lowerMap(), 896 867 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 897 /// \ref supplyType() , \ref flowMap() and \ref potentialMap().868 /// \ref supplyType(). 898 869 /// For example, 899 870 /// \code … … 907 878 /// \ref reset() is called, thus only the modified parameters 908 879 /// have to be set again. See \ref reset() for examples. 880 /// However the underlying digraph must not be modified after this 881 /// class have been constructed, since it copies and extends the graph. 909 882 /// 910 883 /// \param pivot_rule The pivot rule that will be used during the … … 929 902 /// This function resets all the paramaters that have been given 930 903 /// before using functions \ref lowerMap(), \ref upperMap(), 931 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(), 932 /// \ref flowMap() and \ref potentialMap(). 904 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). 933 905 /// 934 906 /// It is useful for multiple run() calls. If this function is not 935 907 /// used, all the parameters given before are kept for the next 936 908 /// \ref run() call. 909 /// However the underlying digraph must not be modified after this 910 /// class have been constructed, since it copies and extends the graph. 937 911 /// 938 912 /// For example, … … 958 932 /// \return <tt>(*this)</tt> 959 933 NetworkSimplex& reset() { 960 delete _plower;961 delete _pupper;962 delete _pcost;963 delete _psupply;964 _plower = NULL;965 _pupper = NULL;966 _pcost = NULL;967 _psupply = NULL;968 _ pstsup= false;934 for (int i = 0; i != _node_num; ++i) { 935 _supply[i] = 0; 936 } 937 for (int i = 0; i != _arc_num; ++i) { 938 _lower[i] = 0; 939 _upper[i] = INF; 940 _cost[i] = 1; 941 } 942 _have_lower = false; 969 943 _stype = GEQ; 970 if (_local_flow) delete _flow_map;971 if (_local_potential) delete _potential_map;972 _flow_map = NULL;973 _potential_map = NULL;974 _local_flow = false;975 _local_potential = false;976 977 944 return *this; 978 945 } … … 1002 969 /// 1003 970 /// \pre \ref run() must be called before using this function. 1004 template <typename Value> 1005 Value totalCost() const { 1006 Value c = 0; 1007 if (_pcost) { 1008 for (ArcIt e(_graph); e != INVALID; ++e) 1009 c += (*_flow_map)[e] * (*_pcost)[e]; 1010 } else { 1011 for (ArcIt e(_graph); e != INVALID; ++e) 1012 c += (*_flow_map)[e]; 971 template <typename Number> 972 Number totalCost() const { 973 Number c = 0; 974 for (ArcIt a(_graph); a != INVALID; ++a) { 975 int i = _arc_id[a]; 976 c += Number(_flow[i]) * Number(_cost[i]); 1013 977 } 1014 978 return c; … … 1027 991 /// \pre \ref run() must be called before using this function. 1028 992 Value flow(const Arc& a) const { 1029 return (*_flow_map)[a]; 1030 } 1031 1032 /// \brief Return a const reference to the flow map. 1033 /// 1034 /// This function returns a const reference to an arc map storing 1035 /// the found flow. 993 return _flow[_arc_id[a]]; 994 } 995 996 /// \brief Return the flow map (the primal solution). 997 /// 998 /// This function copies the flow value on each arc into the given 999 /// map. The \c Value type of the algorithm must be convertible to 1000 /// the \c Value type of the map. 1036 1001 /// 1037 1002 /// \pre \ref run() must be called before using this function. 1038 const FlowMap& flowMap() const { 1039 return *_flow_map; 1003 template <typename FlowMap> 1004 void flowMap(FlowMap &map) const { 1005 for (ArcIt a(_graph); a != INVALID; ++a) { 1006 map.set(a, _flow[_arc_id[a]]); 1007 } 1040 1008 } 1041 1009 … … 1047 1015 /// \pre \ref run() must be called before using this function. 1048 1016 Cost potential(const Node& n) const { 1049 return (*_potential_map)[n];1050 } 1051 1052 /// \brief Return a const reference to the potential map1053 /// (the dual solution).1054 /// 1055 /// This function returns a const reference to a node map storing1056 /// the found potentials, which form the dual solution ofthe1057 /// \ ref min_cost_flow "minimum cost flow problem".1017 return _pi[_node_id[n]]; 1018 } 1019 1020 /// \brief Return the potential map (the dual solution). 1021 /// 1022 /// This function copies the potential (dual value) of each node 1023 /// into the given map. 1024 /// The \c Cost type of the algorithm must be convertible to the 1025 /// \c Value type of the map. 1058 1026 /// 1059 1027 /// \pre \ref run() must be called before using this function. 1060 const PotentialMap& potentialMap() const { 1061 return *_potential_map; 1028 template <typename PotentialMap> 1029 void potentialMap(PotentialMap &map) const { 1030 for (NodeIt n(_graph); n != INVALID; ++n) { 1031 map.set(n, _pi[_node_id[n]]); 1032 } 1062 1033 } 1063 1034 … … 1068 1039 // Initialize internal data structures 1069 1040 bool init() { 1070 // Initialize result maps1071 if (!_flow_map) {1072 _flow_map = new FlowMap(_graph);1073 _local_flow = true;1074 }1075 if (!_potential_map) {1076 _potential_map = new PotentialMap(_graph);1077 _local_potential = true;1078 }1079 1080 // Initialize vectors1081 _node_num = countNodes(_graph);1082 _arc_num = countArcs(_graph);1083 int all_node_num = _node_num + 1;1084 int all_arc_num = _arc_num + _node_num;1085 1041 if (_node_num == 0) return false; 1086 1042 1087 _arc_ref.resize(_arc_num); 1088 _source.resize(all_arc_num); 1089 _target.resize(all_arc_num); 1090 1091 _cap.resize(all_arc_num); 1092 _cost.resize(all_arc_num); 1093 _supply.resize(all_node_num); 1094 _flow.resize(all_arc_num); 1095 _pi.resize(all_node_num); 1096 1097 _parent.resize(all_node_num); 1098 _pred.resize(all_node_num); 1099 _forward.resize(all_node_num); 1100 _thread.resize(all_node_num); 1101 _rev_thread.resize(all_node_num); 1102 _succ_num.resize(all_node_num); 1103 _last_succ.resize(all_node_num); 1104 _state.resize(all_arc_num); 1105 1106 // Initialize node related data 1107 bool valid_supply = true; 1043 // Check the sum of supply values 1108 1044 _sum_supply = 0; 1109 if (!_pstsup && !_psupply) { 1110 _pstsup = true; 1111 _psource = _ptarget = NodeIt(_graph); 1112 _pstflow = 0; 1113 } 1114 if (_psupply) { 1115 int i = 0; 1116 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 1117 _node_id[n] = i; 1118 _supply[i] = (*_psupply)[n]; 1119 _sum_supply += _supply[i]; 1120 } 1121 valid_supply = (_stype == GEQ && _sum_supply <= 0)  1122 (_stype == LEQ && _sum_supply >= 0); 1045 for (int i = 0; i != _node_num; ++i) { 1046 _sum_supply += _supply[i]; 1047 } 1048 if ( !(_stype == GEQ && _sum_supply <= 0  1049 _stype == LEQ && _sum_supply >= 0) ) return false; 1050 1051 // Remove nonzero lower bounds 1052 if (_have_lower) { 1053 for (int i = 0; i != _arc_num; ++i) { 1054 Value c = _lower[i]; 1055 if (c >= 0) { 1056 _cap[i] = _upper[i] < INF ? _upper[i]  c : INF; 1057 } else { 1058 _cap[i] = _upper[i] < INF + c ? _upper[i]  c : INF; 1059 } 1060 _supply[_source[i]] = c; 1061 _supply[_target[i]] += c; 1062 } 1123 1063 } else { 1124 int i = 0; 1125 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 1126 _node_id[n] = i; 1127 _supply[i] = 0; 1128 } 1129 _supply[_node_id[_psource]] = _pstflow; 1130 _supply[_node_id[_ptarget]] = _pstflow; 1131 } 1132 if (!valid_supply) return false; 1064 for (int i = 0; i != _arc_num; ++i) { 1065 _cap[i] = _upper[i]; 1066 } 1067 } 1133 1068 1134 1069 // Initialize artifical cost … … 1144 1079 } 1145 1080 1081 // Initialize arc maps 1082 for (int i = 0; i != _arc_num; ++i) { 1083 _flow[i] = 0; 1084 _state[i] = STATE_LOWER; 1085 } 1086 1146 1087 // Set data for the artificial root node 1147 1088 _root = _node_num; … … 1150 1091 _thread[_root] = 0; 1151 1092 _rev_thread[0] = _root; 1152 _succ_num[_root] = all_node_num;1093 _succ_num[_root] = _node_num + 1; 1153 1094 _last_succ[_root] = _root  1; 1154 1095 _supply[_root] = _sum_supply; 1155 if (_sum_supply < 0) { 1156 _pi[_root] = ART_COST; 1157 } else { 1158 _pi[_root] = ART_COST; 1159 } 1160 1161 // Store the arcs in a mixed order 1162 int k = std::max(int(std::sqrt(double(_arc_num))), 10); 1163 int i = 0; 1164 for (ArcIt e(_graph); e != INVALID; ++e) { 1165 _arc_ref[i] = e; 1166 if ((i += k) >= _arc_num) i = (i % k) + 1; 1167 } 1168 1169 // Initialize arc maps 1170 if (_pupper && _pcost) { 1171 for (int i = 0; i != _arc_num; ++i) { 1172 Arc e = _arc_ref[i]; 1173 _source[i] = _node_id[_graph.source(e)]; 1174 _target[i] = _node_id[_graph.target(e)]; 1175 _cap[i] = (*_pupper)[e]; 1176 _cost[i] = (*_pcost)[e]; 1177 _flow[i] = 0; 1178 _state[i] = STATE_LOWER; 1179 } 1180 } else { 1181 for (int i = 0; i != _arc_num; ++i) { 1182 Arc e = _arc_ref[i]; 1183 _source[i] = _node_id[_graph.source(e)]; 1184 _target[i] = _node_id[_graph.target(e)]; 1185 _flow[i] = 0; 1186 _state[i] = STATE_LOWER; 1187 } 1188 if (_pupper) { 1189 for (int i = 0; i != _arc_num; ++i) 1190 _cap[i] = (*_pupper)[_arc_ref[i]]; 1191 } else { 1192 for (int i = 0; i != _arc_num; ++i) 1193 _cap[i] = INF; 1194 } 1195 if (_pcost) { 1196 for (int i = 0; i != _arc_num; ++i) 1197 _cost[i] = (*_pcost)[_arc_ref[i]]; 1198 } else { 1199 for (int i = 0; i != _arc_num; ++i) 1200 _cost[i] = 1; 1201 } 1202 } 1203 1204 // Remove nonzero lower bounds 1205 if (_plower) { 1206 for (int i = 0; i != _arc_num; ++i) { 1207 Value c = (*_plower)[_arc_ref[i]]; 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 } 1219 _supply[_source[i]] = c; 1220 _supply[_target[i]] += c; 1221 } 1222 } 1223 } 1096 _pi[_root] = _sum_supply < 0 ? ART_COST : ART_COST; 1224 1097 1225 1098 // Add artificial arcs and initialize the spanning tree data structure 1226 1099 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1100 _parent[u] = _root; 1101 _pred[u] = e; 1227 1102 _thread[u] = u + 1; 1228 1103 _rev_thread[u + 1] = u; 1229 1104 _succ_num[u] = 1; 1230 1105 _last_succ[u] = u; 1231 _parent[u] = _root;1232 _pred[u] = e;1233 1106 _cost[e] = ART_COST; 1234 1107 _cap[e] = INF; … … 1518 1391 } 1519 1392 1520 // Copy flow values to _flow_map1521 if (_ plower) {1393 // Transform the solution and the supply map to the original form 1394 if (_have_lower) { 1522 1395 for (int i = 0; i != _arc_num; ++i) { 1523 Arc e = _arc_ref[i]; 1524 _flow_map>set(e, (*_plower)[e] + _flow[i]); 1525 } 1526 } else { 1527 for (int i = 0; i != _arc_num; ++i) { 1528 _flow_map>set(_arc_ref[i], _flow[i]); 1529 } 1530 } 1531 // Copy potential values to _potential_map 1532 for (NodeIt n(_graph); n != INVALID; ++n) { 1533 _potential_map>set(n, _pi[_node_id[n]]); 1396 Value c = _lower[i]; 1397 if (c != 0) { 1398 _flow[i] += c; 1399 _supply[_source[i]] += c; 1400 _supply[_target[i]] = c; 1401 } 1402 } 1534 1403 } 1535 1404
Note: See TracChangeset
for help on using the changeset viewer.