0
2
0
149
280
... | ... |
@@ -72,21 +72,10 @@ |
72 | 72 |
{ |
73 | 73 |
public: |
74 | 74 |
|
75 |
/// The |
|
75 |
/// The type of the flow amounts, capacity bounds and supply values |
|
76 | 76 |
typedef V Value; |
77 |
/// The |
|
77 |
/// The type of the arc costs |
|
78 | 78 |
typedef C Cost; |
79 |
#ifdef DOXYGEN |
|
80 |
/// The type of the flow map |
|
81 |
typedef GR::ArcMap<Value> FlowMap; |
|
82 |
/// The type of the potential map |
|
83 |
typedef GR::NodeMap<Cost> PotentialMap; |
|
84 |
#else |
|
85 |
/// The type of the flow map |
|
86 |
typedef typename GR::template ArcMap<Value> FlowMap; |
|
87 |
/// The type of the potential map |
|
88 |
typedef typename GR::template NodeMap<Cost> PotentialMap; |
|
89 |
#endif |
|
90 | 79 |
|
91 | 80 |
public: |
92 | 81 |
|
... | ... |
@@ -206,15 +195,11 @@ |
206 | 195 |
|
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> |
|
202 |
typedef std::vector<Value> ValueVector; |
|
218 | 203 |
typedef std::vector<Cost> CostVector; |
219 | 204 |
|
220 | 205 |
// State constants for arcs |
... | ... |
@@ -232,34 +217,23 @@ |
232 | 217 |
int _arc_num; |
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 |
|
|
220 |
bool _have_lower; |
|
242 | 221 |
SupplyType _stype; |
243 |
|
|
244 | 222 |
Value _sum_supply; |
245 | 223 |
|
246 |
// Result maps |
|
247 |
FlowMap *_flow_map; |
|
248 |
PotentialMap *_potential_map; |
|
249 |
bool _local_flow; |
|
250 |
bool _local_potential; |
|
251 |
|
|
252 | 224 |
// Data structures for storing the digraph |
253 | 225 |
IntNodeMap _node_id; |
254 |
|
|
226 |
IntArcMap _arc_id; |
|
255 | 227 |
IntVector _source; |
256 | 228 |
IntVector _target; |
257 | 229 |
|
258 | 230 |
// Node and arc data |
259 |
|
|
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 |
|
265 | 239 |
// Data for storing the spanning tree structure |
... | ... |
@@ -689,12 +663,7 @@ |
689 | 663 |
/// |
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() : |
700 | 669 |
std::numeric_limits<Value>::max()) |
... | ... |
@@ -704,12 +673,58 @@ |
704 | 673 |
"The flow type of NetworkSimplex must be signed"); |
705 | 674 |
LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
706 | 675 |
"The cost type of NetworkSimplex must be signed"); |
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; |
|
707 | 715 |
} |
708 | 716 |
|
709 |
/// Destructor. |
|
710 |
~NetworkSimplex() { |
|
711 |
if (_local_flow) delete _flow_map; |
|
712 |
if (_local_potential) delete _potential_map; |
|
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 | 730 |
/// \name Parameters |
... | ... |
@@ -731,10 +746,9 @@ |
731 | 746 |
/// \return <tt>(*this)</tt> |
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 |
|
|
751 |
_lower[_arc_id[a]] = map[a]; |
|
738 | 752 |
} |
739 | 753 |
return *this; |
740 | 754 |
} |
... | ... |
@@ -753,10 +767,8 @@ |
753 | 767 |
/// \return <tt>(*this)</tt> |
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 |
|
|
771 |
_upper[_arc_id[a]] = map[a]; |
|
760 | 772 |
} |
761 | 773 |
return *this; |
762 | 774 |
} |
... | ... |
@@ -774,10 +786,8 @@ |
774 | 786 |
/// \return <tt>(*this)</tt> |
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 |
|
|
790 |
_cost[_arc_id[a]] = map[a]; |
|
781 | 791 |
} |
782 | 792 |
return *this; |
783 | 793 |
} |
... | ... |
@@ -796,11 +806,8 @@ |
796 | 806 |
/// \return <tt>(*this)</tt> |
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 |
|
|
810 |
_supply[_node_id[n]] = map[n]; |
|
804 | 811 |
} |
805 | 812 |
return *this; |
806 | 813 |
} |
... | ... |
@@ -824,12 +831,11 @@ |
824 | 831 |
/// |
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 |
} |
835 | 841 |
|
... | ... |
@@ -847,41 +853,6 @@ |
847 | 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 | 858 |
/// \name Execution Control |
... | ... |
@@ -894,7 +865,7 @@ |
894 | 865 |
/// This function runs the algorithm. |
895 | 866 |
/// The paramters can be specified using functions \ref lowerMap(), |
896 | 867 |
/// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), |
897 |
/// \ref supplyType() |
|
868 |
/// \ref supplyType(). |
|
898 | 869 |
/// For example, |
899 | 870 |
/// \code |
900 | 871 |
/// NetworkSimplex<ListDigraph> ns(graph); |
... | ... |
@@ -906,6 +877,8 @@ |
906 | 877 |
/// that have been given are kept for the next call, unless |
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 |
911 | 884 |
/// algorithm. For more information see \ref PivotRule. |
... | ... |
@@ -928,12 +901,13 @@ |
928 | 901 |
/// |
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, |
939 | 913 |
/// \code |
... | ... |
@@ -957,23 +931,16 @@ |
957 | 931 |
/// |
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 |
|
|
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 |
} |
979 | 946 |
|
... | ... |
@@ -1001,15 +968,12 @@ |
1001 | 968 |
/// function. |
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 |
|
|
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; |
1015 | 979 |
} |
... | ... |
@@ -1026,17 +990,21 @@ |
1026 | 990 |
/// |
1027 | 991 |
/// \pre \ref run() must be called before using this function. |
1028 | 992 |
Value flow(const Arc& a) const { |
1029 |
return |
|
993 |
return _flow[_arc_id[a]]; |
|
1030 | 994 |
} |
1031 | 995 |
|
1032 |
/// \brief Return |
|
996 |
/// \brief Return the flow map (the primal solution). |
|
1033 | 997 |
/// |
1034 |
/// This function returns a const reference to an arc map storing |
|
1035 |
/// the found flow. |
|
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 |
|
1042 | 1010 |
/// \brief Return the potential (dual value) of the given node. |
... | ... |
@@ -1046,19 +1014,22 @@ |
1046 | 1014 |
/// |
1047 | 1015 |
/// \pre \ref run() must be called before using this function. |
1048 | 1016 |
Cost potential(const Node& n) const { |
1049 |
return |
|
1017 |
return _pi[_node_id[n]]; |
|
1050 | 1018 |
} |
1051 | 1019 |
|
1052 |
/// \brief Return a const reference to the potential map |
|
1053 |
/// (the dual solution). |
|
1020 |
/// \brief Return the potential map (the dual solution). |
|
1054 | 1021 |
/// |
1055 |
/// This function returns a const reference to a node map storing |
|
1056 |
/// the found potentials, which form the dual solution of the |
|
1057 |
/// |
|
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 |
|
1064 | 1035 |
/// @} |
... | ... |
@@ -1067,69 +1038,33 @@ |
1067 | 1038 |
|
1068 | 1039 |
// Initialize internal data structures |
1069 | 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 | 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 |
|
|
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]; |
|
1045 |
for (int i = 0; i != _node_num; ++i) { |
|
1119 | 1046 |
_sum_supply += _supply[i]; |
1120 | 1047 |
} |
1121 |
valid_supply = (_stype == GEQ && _sum_supply <= 0) || |
|
1122 |
(_stype == LEQ && _sum_supply >= 0); |
|
1048 |
if ( !(_stype == GEQ && _sum_supply <= 0 || |
|
1049 |
_stype == LEQ && _sum_supply >= 0) ) return false; |
|
1050 |
|
|
1051 |
// Remove non-zero 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; |
|
1123 | 1057 |
} else { |
1124 |
int i = 0; |
|
1125 |
for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
|
1126 |
_node_id[n] = i; |
|
1127 |
_supply[i] = 0; |
|
1058 |
_cap[i] = _upper[i] < INF + c ? _upper[i] - c : INF; |
|
1128 | 1059 |
} |
1129 |
_supply[_node_id[_psource]] = _pstflow; |
|
1130 |
_supply[_node_id[_ptarget]] = -_pstflow; |
|
1060 |
_supply[_source[i]] -= c; |
|
1061 |
_supply[_target[i]] += c; |
|
1131 | 1062 |
} |
1132 |
|
|
1063 |
} else { |
|
1064 |
for (int i = 0; i != _arc_num; ++i) { |
|
1065 |
_cap[i] = _upper[i]; |
|
1066 |
} |
|
1067 |
} |
|
1133 | 1068 |
|
1134 | 1069 |
// Initialize artifical cost |
1135 | 1070 |
Cost ART_COST; |
... | ... |
@@ -1143,93 +1078,31 @@ |
1143 | 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 | 1087 |
// Set data for the artificial root node |
1147 | 1088 |
_root = _node_num; |
1148 | 1089 |
_parent[_root] = -1; |
1149 | 1090 |
_pred[_root] = -1; |
1150 | 1091 |
_thread[_root] = 0; |
1151 | 1092 |
_rev_thread[0] = _root; |
1152 |
_succ_num[_root] = |
|
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 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 |
|
|
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; |
1235 | 1108 |
_state[e] = STATE_TREE; |
... | ... |
@@ -1517,21 +1390,17 @@ |
1517 | 1390 |
} |
1518 | 1391 |
} |
1519 | 1392 |
|
1520 |
// Copy flow values to _flow_map |
|
1521 |
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]); |
|
1396 |
Value c = _lower[i]; |
|
1397 |
if (c != 0) { |
|
1398 |
_flow[i] += c; |
|
1399 |
_supply[_source[i]] += c; |
|
1400 |
_supply[_target[i]] -= c; |
|
1525 | 1401 |
} |
1526 |
} else { |
|
1527 |
for (int i = 0; i != _arc_num; ++i) { |
|
1528 |
_flow_map->set(_arc_ref[i], _flow[i]); |
|
1529 | 1402 |
} |
1530 | 1403 |
} |
1531 |
// Copy potential values to _potential_map |
|
1532 |
for (NodeIt n(_graph); n != INVALID; ++n) { |
|
1533 |
_potential_map->set(n, _pi[_node_id[n]]); |
|
1534 |
} |
|
1535 | 1404 |
|
1536 | 1405 |
return OPTIMAL; |
1537 | 1406 |
} |
... | ... |
@@ -84,7 +84,7 @@ |
84 | 84 |
}; |
85 | 85 |
|
86 | 86 |
// Check the interface of an MCF algorithm |
87 |
template <typename GR, typename |
|
87 |
template <typename GR, typename Value, typename Cost> |
|
88 | 88 |
class McfClassConcept |
89 | 89 |
{ |
90 | 90 |
public: |
... | ... |
@@ -95,6 +95,7 @@ |
95 | 95 |
checkConcept<concepts::Digraph, GR>(); |
96 | 96 |
|
97 | 97 |
MCF mcf(g); |
98 |
const MCF& const_mcf = mcf; |
|
98 | 99 |
|
99 | 100 |
b = mcf.reset() |
100 | 101 |
.lowerMap(lower) |
... | ... |
@@ -102,47 +103,38 @@ |
102 | 103 |
.costMap(cost) |
103 | 104 |
.supplyMap(sup) |
104 | 105 |
.stSupply(n, n, k) |
105 |
.flowMap(flow) |
|
106 |
.potentialMap(pot) |
|
107 | 106 |
.run(); |
108 | 107 |
|
109 |
const MCF& const_mcf = mcf; |
|
110 |
|
|
111 |
const typename MCF::FlowMap &fm = const_mcf.flowMap(); |
|
112 |
const typename MCF::PotentialMap &pm = const_mcf.potentialMap(); |
|
113 |
|
|
114 | 108 |
c = const_mcf.totalCost(); |
115 |
|
|
109 |
x = const_mcf.template totalCost<double>(); |
|
116 | 110 |
v = const_mcf.flow(a); |
117 | 111 |
c = const_mcf.potential(n); |
118 |
|
|
119 |
v = const_mcf.INF; |
|
120 |
|
|
121 |
ignore_unused_variable_warning(fm); |
|
122 |
ignore_unused_variable_warning(pm); |
|
123 |
ignore_unused_variable_warning(x); |
|
112 |
const_mcf.flowMap(fm); |
|
113 |
const_mcf.potentialMap(pm); |
|
124 | 114 |
} |
125 | 115 |
|
126 | 116 |
typedef typename GR::Node Node; |
127 | 117 |
typedef typename GR::Arc Arc; |
128 |
typedef concepts::ReadMap<Node, Flow> NM; |
|
129 |
typedef concepts::ReadMap<Arc, Flow> FAM; |
|
118 |
typedef concepts::ReadMap<Node, Value> NM; |
|
119 |
typedef concepts::ReadMap<Arc, Value> VAM; |
|
130 | 120 |
typedef concepts::ReadMap<Arc, Cost> CAM; |
121 |
typedef concepts::WriteMap<Arc, Value> FlowMap; |
|
122 |
typedef concepts::WriteMap<Node, Cost> PotMap; |
|
131 | 123 |
|
132 | 124 |
const GR &g; |
133 |
const FAM &lower; |
|
134 |
const FAM &upper; |
|
125 |
const VAM &lower; |
|
126 |
const VAM &upper; |
|
135 | 127 |
const CAM &cost; |
136 | 128 |
const NM ⊃ |
137 | 129 |
const Node &n; |
138 | 130 |
const Arc &a; |
139 |
const Flow &k; |
|
140 |
Flow v; |
|
141 |
|
|
131 |
const Value &k; |
|
132 |
FlowMap fm; |
|
133 |
PotMap pm; |
|
142 | 134 |
bool b; |
143 |
|
|
144 |
typename MCF::FlowMap &flow; |
|
145 |
|
|
135 |
double x; |
|
136 |
typename MCF::Value v; |
|
137 |
typename MCF::Cost c; |
|
146 | 138 |
}; |
147 | 139 |
|
148 | 140 |
}; |
... | ... |
@@ -221,11 +213,14 @@ |
221 | 213 |
{ |
222 | 214 |
check(mcf_result == result, "Wrong result " + test_id); |
223 | 215 |
if (optimal) { |
224 |
|
|
216 |
typename GR::template ArcMap<typename SM::Value> flow(gr); |
|
217 |
typename GR::template NodeMap<typename CM::Value> pi(gr); |
|
218 |
mcf.flowMap(flow); |
|
219 |
mcf.potentialMap(pi); |
|
220 |
check(checkFlow(gr, lower, upper, supply, flow, type), |
|
225 | 221 |
"The flow is not feasible " + test_id); |
226 | 222 |
check(mcf.totalCost() == total, "The flow is not optimal " + test_id); |
227 |
check(checkPotential(gr, lower, upper, cost, supply, mcf.flowMap(), |
|
228 |
mcf.potentialMap()), |
|
223 |
check(checkPotential(gr, lower, upper, cost, supply, flow, pi), |
|
229 | 224 |
"Wrong potentials " + test_id); |
230 | 225 |
} |
231 | 226 |
} |
... | ... |
@@ -234,11 +229,13 @@ |
234 | 229 |
{ |
235 | 230 |
// Check the interfaces |
236 | 231 |
{ |
237 |
typedef int Flow; |
|
238 |
typedef int Cost; |
|
239 | 232 |
typedef concepts::Digraph GR; |
240 |
checkConcept< McfClassConcept<GR, Flow, Cost>, |
|
241 |
NetworkSimplex<GR, Flow, Cost> >(); |
|
233 |
checkConcept< McfClassConcept<GR, int, int>, |
|
234 |
NetworkSimplex<GR> >(); |
|
235 |
checkConcept< McfClassConcept<GR, double, double>, |
|
236 |
NetworkSimplex<GR, double> >(); |
|
237 |
checkConcept< McfClassConcept<GR, int, double>, |
|
238 |
NetworkSimplex<GR, int, double> >(); |
|
242 | 239 |
} |
243 | 240 |
|
244 | 241 |
// Run various MCF tests |
0 comments (0 inline)