687 /// |
661 /// |
688 /// The constructor of the class. |
662 /// The constructor of the class. |
689 /// |
663 /// |
690 /// \param graph The digraph the algorithm runs on. |
664 /// \param graph The digraph the algorithm runs on. |
691 NetworkSimplex(const GR& graph) : |
665 NetworkSimplex(const GR& graph) : |
692 _graph(graph), |
666 _graph(graph), _node_id(graph), _arc_id(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), |
|
698 INF(std::numeric_limits<Value>::has_infinity ? |
667 INF(std::numeric_limits<Value>::has_infinity ? |
699 std::numeric_limits<Value>::infinity() : |
668 std::numeric_limits<Value>::infinity() : |
700 std::numeric_limits<Value>::max()) |
669 std::numeric_limits<Value>::max()) |
701 { |
670 { |
702 // Check the value types |
671 // Check the value types |
703 LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
672 LEMON_ASSERT(std::numeric_limits<Value>::is_signed, |
704 "The flow type of NetworkSimplex must be signed"); |
673 "The flow type of NetworkSimplex must be signed"); |
705 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
674 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
706 "The cost type of NetworkSimplex must be signed"); |
675 "The cost type of NetworkSimplex must be signed"); |
707 } |
676 |
708 |
677 // Resize vectors |
709 /// Destructor. |
678 _node_num = countNodes(_graph); |
710 ~NetworkSimplex() { |
679 _arc_num = countArcs(_graph); |
711 if (_local_flow) delete _flow_map; |
680 int all_node_num = _node_num + 1; |
712 if (_local_potential) delete _potential_map; |
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 |
715 /// \name Parameters |
730 /// \name Parameters |
716 /// The parameters of the algorithm can be specified using these |
731 /// The parameters of the algorithm can be specified using these |
717 /// functions. |
732 /// functions. |
845 NetworkSimplex& supplyType(SupplyType supply_type) { |
851 NetworkSimplex& supplyType(SupplyType supply_type) { |
846 _stype = supply_type; |
852 _stype = supply_type; |
847 return *this; |
853 return *this; |
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 will |
|
854 /// be allocated automatically. The destructor deallocates this |
|
855 /// 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 storing |
|
870 /// the dual solution. |
|
871 /// If it is not used before calling \ref run(), an instance will |
|
872 /// be allocated automatically. The destructor deallocates this |
|
873 /// 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 |
887 /// \name Execution Control |
858 /// \name Execution Control |
888 /// The algorithm can be executed using \ref run(). |
859 /// The algorithm can be executed using \ref run(). |
889 |
860 |
926 |
899 |
927 /// \brief Reset all the parameters that have been given before. |
900 /// \brief Reset all the parameters that have been given before. |
928 /// |
901 /// |
929 /// This function resets all the paramaters that have been given |
902 /// This function resets all the paramaters that have been given |
930 /// before using functions \ref lowerMap(), \ref upperMap(), |
903 /// before using functions \ref lowerMap(), \ref upperMap(), |
931 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(), |
904 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). |
932 /// \ref flowMap() and \ref potentialMap(). |
|
933 /// |
905 /// |
934 /// It is useful for multiple run() calls. If this function is not |
906 /// It is useful for multiple run() calls. If this function is not |
935 /// used, all the parameters given before are kept for the next |
907 /// used, all the parameters given before are kept for the next |
936 /// \ref run() call. |
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 /// For example, |
912 /// For example, |
939 /// \code |
913 /// \code |
940 /// NetworkSimplex<ListDigraph> ns(graph); |
914 /// NetworkSimplex<ListDigraph> ns(graph); |
941 /// |
915 /// |
1024 /// |
988 /// |
1025 /// This function returns the flow on the given arc. |
989 /// This function returns the flow on the given arc. |
1026 /// |
990 /// |
1027 /// \pre \ref run() must be called before using this function. |
991 /// \pre \ref run() must be called before using this function. |
1028 Value flow(const Arc& a) const { |
992 Value flow(const Arc& a) const { |
1029 return (*_flow_map)[a]; |
993 return _flow[_arc_id[a]]; |
1030 } |
994 } |
1031 |
995 |
1032 /// \brief Return a const reference to the flow map. |
996 /// \brief Return the flow map (the primal solution). |
1033 /// |
997 /// |
1034 /// This function returns a const reference to an arc map storing |
998 /// This function copies the flow value on each arc into the given |
1035 /// the found flow. |
999 /// map. The \c Value type of the algorithm must be convertible to |
|
1000 /// the \c Value type of the map. |
1036 /// |
1001 /// |
1037 /// \pre \ref run() must be called before using this function. |
1002 /// \pre \ref run() must be called before using this function. |
1038 const FlowMap& flowMap() const { |
1003 template <typename FlowMap> |
1039 return *_flow_map; |
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 |
1042 /// \brief Return the potential (dual value) of the given node. |
1010 /// \brief Return the potential (dual value) of the given node. |
1043 /// |
1011 /// |
1044 /// This function returns the potential (dual value) of the |
1012 /// This function returns the potential (dual value) of the |
1045 /// given node. |
1013 /// given node. |
1046 /// |
1014 /// |
1047 /// \pre \ref run() must be called before using this function. |
1015 /// \pre \ref run() must be called before using this function. |
1048 Cost potential(const Node& n) const { |
1016 Cost potential(const Node& n) const { |
1049 return (*_potential_map)[n]; |
1017 return _pi[_node_id[n]]; |
1050 } |
1018 } |
1051 |
1019 |
1052 /// \brief Return a const reference to the potential map |
1020 /// \brief Return the potential map (the dual solution). |
1053 /// (the dual solution). |
1021 /// |
1054 /// |
1022 /// This function copies the potential (dual value) of each node |
1055 /// This function returns a const reference to a node map storing |
1023 /// into the given map. |
1056 /// the found potentials, which form the dual solution of the |
1024 /// The \c Cost type of the algorithm must be convertible to the |
1057 /// \ref min_cost_flow "minimum cost flow problem". |
1025 /// \c Value type of the map. |
1058 /// |
1026 /// |
1059 /// \pre \ref run() must be called before using this function. |
1027 /// \pre \ref run() must be called before using this function. |
1060 const PotentialMap& potentialMap() const { |
1028 template <typename PotentialMap> |
1061 return *_potential_map; |
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 |
1064 /// @} |
1035 /// @} |
1065 |
1036 |
1066 private: |
1037 private: |
1067 |
1038 |
1068 // Initialize internal data structures |
1039 // Initialize internal data structures |
1069 bool init() { |
1040 bool init() { |
1070 // Initialize result maps |
|
1071 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 vectors |
|
1081 _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 if (_node_num == 0) return false; |
1041 if (_node_num == 0) return false; |
1086 |
1042 |
1087 _arc_ref.resize(_arc_num); |
1043 // Check the sum of supply values |
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; |
|
1108 _sum_supply = 0; |
1044 _sum_supply = 0; |
1109 if (!_pstsup && !_psupply) { |
1045 for (int i = 0; i != _node_num; ++i) { |
1110 _pstsup = true; |
1046 _sum_supply += _supply[i]; |
1111 _psource = _ptarget = NodeIt(_graph); |
1047 } |
1112 _pstflow = 0; |
1048 if ( !(_stype == GEQ && _sum_supply <= 0 || |
1113 } |
1049 _stype == LEQ && _sum_supply >= 0) ) return false; |
1114 if (_psupply) { |
1050 |
1115 int i = 0; |
1051 // Remove non-zero lower bounds |
1116 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
1052 if (_have_lower) { |
1117 _node_id[n] = i; |
1053 for (int i = 0; i != _arc_num; ++i) { |
1118 _supply[i] = (*_psupply)[n]; |
1054 Value c = _lower[i]; |
1119 _sum_supply += _supply[i]; |
1055 if (c >= 0) { |
1120 } |
1056 _cap[i] = _upper[i] < INF ? _upper[i] - c : INF; |
1121 valid_supply = (_stype == GEQ && _sum_supply <= 0) || |
1057 } else { |
1122 (_stype == LEQ && _sum_supply >= 0); |
1058 _cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF; |
|
1059 } |
|
1060 _supply[_source[i]] -= c; |
|
1061 _supply[_target[i]] += c; |
|
1062 } |
1123 } else { |
1063 } else { |
1124 int i = 0; |
1064 for (int i = 0; i != _arc_num; ++i) { |
1125 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
1065 _cap[i] = _upper[i]; |
1126 _node_id[n] = i; |
1066 } |
1127 _supply[i] = 0; |
1067 } |
1128 } |
|
1129 _supply[_node_id[_psource]] = _pstflow; |
|
1130 _supply[_node_id[_ptarget]] = -_pstflow; |
|
1131 } |
|
1132 if (!valid_supply) return false; |
|
1133 |
1068 |
1134 // Initialize artifical cost |
1069 // Initialize artifical cost |
1135 Cost ART_COST; |
1070 Cost ART_COST; |
1136 if (std::numeric_limits<Cost>::is_exact) { |
1071 if (std::numeric_limits<Cost>::is_exact) { |
1137 ART_COST = std::numeric_limits<Cost>::max() / 4 + 1; |
1072 ART_COST = std::numeric_limits<Cost>::max() / 4 + 1; |
1141 if (_cost[i] > ART_COST) ART_COST = _cost[i]; |
1076 if (_cost[i] > ART_COST) ART_COST = _cost[i]; |
1142 } |
1077 } |
1143 ART_COST = (ART_COST + 1) * _node_num; |
1078 ART_COST = (ART_COST + 1) * _node_num; |
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 // Set data for the artificial root node |
1087 // Set data for the artificial root node |
1147 _root = _node_num; |
1088 _root = _node_num; |
1148 _parent[_root] = -1; |
1089 _parent[_root] = -1; |
1149 _pred[_root] = -1; |
1090 _pred[_root] = -1; |
1150 _thread[_root] = 0; |
1091 _thread[_root] = 0; |
1151 _rev_thread[0] = _root; |
1092 _rev_thread[0] = _root; |
1152 _succ_num[_root] = all_node_num; |
1093 _succ_num[_root] = _node_num + 1; |
1153 _last_succ[_root] = _root - 1; |
1094 _last_succ[_root] = _root - 1; |
1154 _supply[_root] = -_sum_supply; |
1095 _supply[_root] = -_sum_supply; |
1155 if (_sum_supply < 0) { |
1096 _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST; |
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 non-zero 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 } |
|
1224 |
1097 |
1225 // Add artificial arcs and initialize the spanning tree data structure |
1098 // Add artificial arcs and initialize the spanning tree data structure |
1226 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
1099 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { |
|
1100 _parent[u] = _root; |
|
1101 _pred[u] = e; |
1227 _thread[u] = u + 1; |
1102 _thread[u] = u + 1; |
1228 _rev_thread[u + 1] = u; |
1103 _rev_thread[u + 1] = u; |
1229 _succ_num[u] = 1; |
1104 _succ_num[u] = 1; |
1230 _last_succ[u] = u; |
1105 _last_succ[u] = u; |
1231 _parent[u] = _root; |
|
1232 _pred[u] = e; |
|
1233 _cost[e] = ART_COST; |
1106 _cost[e] = ART_COST; |
1234 _cap[e] = INF; |
1107 _cap[e] = INF; |
1235 _state[e] = STATE_TREE; |
1108 _state[e] = STATE_TREE; |
1236 if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) { |
1109 if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) { |
1237 _flow[e] = _supply[u]; |
1110 _flow[e] = _supply[u]; |