48 /// simplex method directly for the minimum cost flow problem. |
45 /// simplex method directly for the minimum cost flow problem. |
49 /// It is one of the most efficient solution methods. |
46 /// It is one of the most efficient solution methods. |
50 /// |
47 /// |
51 /// In general this class is the fastest implementation available |
48 /// In general this class is the fastest implementation available |
52 /// in LEMON for the minimum cost flow problem. |
49 /// in LEMON for the minimum cost flow problem. |
53 /// Moreover it supports both direction of the supply/demand inequality |
50 /// Moreover it supports both directions of the supply/demand inequality |
54 /// constraints. For more information see \ref ProblemType. |
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 /// \tparam GR The digraph type the algorithm runs on. |
58 /// \tparam GR The digraph type the algorithm runs on. |
57 /// \tparam F The value type used for flow amounts, capacity bounds |
59 /// \tparam F The value type used for flow amounts, capacity bounds |
58 /// and supply values in the algorithm. By default it is \c int. |
60 /// and supply values in the algorithm. By default it is \c int. |
59 /// \tparam C The value type used for costs and potentials in the |
61 /// \tparam C The value type used for costs and potentials in the |
86 typedef typename GR::template NodeMap<Cost> PotentialMap; |
88 typedef typename GR::template NodeMap<Cost> PotentialMap; |
87 #endif |
89 #endif |
88 |
90 |
89 public: |
91 public: |
90 |
92 |
91 /// \brief Enum type for selecting the pivot rule. |
93 /// \brief Problem type constants for the \c run() function. |
92 /// |
94 /// |
93 /// Enum type for selecting the pivot rule for the \ref run() |
95 /// Enum type containing the problem type constants that can be |
94 /// function. |
96 /// returned by the \ref run() function of the algorithm. |
95 /// |
97 enum ProblemType { |
96 /// \ref NetworkSimplex provides five different pivot rule |
98 /// The problem has no feasible solution (flow). |
97 /// implementations that significantly affect the running time |
99 INFEASIBLE, |
98 /// of the algorithm. |
100 /// The problem has optimal solution (i.e. it is feasible and |
99 /// By default \ref BLOCK_SEARCH "Block Search" is used, which |
101 /// bounded), and the algorithm has found optimal flow and node |
100 /// proved to be the most efficient and the most robust on various |
102 /// potentials (primal and dual solutions). |
101 /// test inputs according to our benchmark tests. |
103 OPTIMAL, |
102 /// However another pivot rule can be selected using the \ref run() |
104 /// The objective function of the problem is unbounded, i.e. |
103 /// function with the proper parameter. |
105 /// there is a directed cycle having negative total cost and |
104 enum PivotRule { |
106 /// infinite upper bound. |
105 |
107 UNBOUNDED |
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 |
|
132 }; |
108 }; |
133 |
109 |
134 /// \brief Enum type for selecting the problem type. |
110 /// \brief Constants for selecting the type of the supply constraints. |
135 /// |
111 /// |
136 /// Enum type for selecting the problem type, i.e. the direction of |
112 /// Enum type containing constants for selecting the supply type, |
137 /// the inequalities in the supply/demand constraints of the |
113 /// i.e. the direction of the inequalities in the supply/demand |
138 /// \ref min_cost_flow "minimum cost flow problem". |
114 /// constraints of the \ref min_cost_flow "minimum cost flow problem". |
139 /// |
115 /// |
140 /// The default problem type is \c GEQ, since this form is supported |
116 /// The default supply type is \c GEQ, since this form is supported |
141 /// by other minimum cost flow algorithms and the \ref Circulation |
117 /// by other minimum cost flow algorithms and the \ref Circulation |
142 /// algorithm as well. |
118 /// algorithm, as well. |
143 /// The \c LEQ problem type can be selected using the \ref problemType() |
119 /// The \c LEQ problem type can be selected using the \ref supplyType() |
144 /// function. |
120 /// function. |
145 /// |
121 /// |
146 /// Note that the equality form is a special case of both problem type. |
122 /// Note that the equality form is a special case of both supply types. |
147 enum ProblemType { |
123 enum SupplyType { |
148 |
124 |
149 /// This option means that there are "<em>greater or equal</em>" |
125 /// This option means that there are <em>"greater or equal"</em> |
150 /// constraints in the defintion, i.e. the exact formulation of the |
126 /// supply/demand constraints in the definition, i.e. the exact |
151 /// problem is the following. |
127 /// formulation of the problem is the following. |
152 /** |
128 /** |
153 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
129 \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 |
130 \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq |
155 sup(u) \quad \forall u\in V \f] |
131 sup(u) \quad \forall u\in V \f] |
156 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
132 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
162 /// satisfied. |
138 /// satisfied. |
163 GEQ, |
139 GEQ, |
164 /// It is just an alias for the \c GEQ option. |
140 /// It is just an alias for the \c GEQ option. |
165 CARRY_SUPPLIES = GEQ, |
141 CARRY_SUPPLIES = GEQ, |
166 |
142 |
167 /// This option means that there are "<em>less or equal</em>" |
143 /// This option means that there are <em>"less or equal"</em> |
168 /// constraints in the defintion, i.e. the exact formulation of the |
144 /// supply/demand constraints in the definition, i.e. the exact |
169 /// problem is the following. |
145 /// formulation of the problem is the following. |
170 /** |
146 /** |
171 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] |
147 \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 |
148 \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq |
173 sup(u) \quad \forall u\in V \f] |
149 sup(u) \quad \forall u\in V \f] |
174 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
150 \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f] |
180 /// nodes. |
156 /// nodes. |
181 LEQ, |
157 LEQ, |
182 /// It is just an alias for the \c LEQ option. |
158 /// It is just an alias for the \c LEQ option. |
183 SATISFY_DEMANDS = LEQ |
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 private: |
205 private: |
187 |
206 |
188 TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
207 TEMPLATE_DIGRAPH_TYPEDEFS(GR); |
189 |
208 |
190 typedef typename GR::template ArcMap<Flow> FlowArcMap; |
209 typedef typename GR::template ArcMap<Flow> FlowArcMap; |
256 // Temporary data used in the current pivot iteration |
277 // Temporary data used in the current pivot iteration |
257 int in_arc, join, u_in, v_in, u_out, v_out; |
278 int in_arc, join, u_in, v_in, u_out, v_out; |
258 int first, second, right, last; |
279 int first, second, right, last; |
259 int stem, par_stem, new_stem; |
280 int stem, par_stem, new_stem; |
260 Flow delta; |
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 private: |
292 private: |
263 |
293 |
264 // Implementation of the First Eligible pivot rule |
294 // Implementation of the First Eligible pivot rule |
265 class FirstEligiblePivotRule |
295 class FirstEligiblePivotRule |
659 /// |
689 /// |
660 /// \param graph The digraph the algorithm runs on. |
690 /// \param graph The digraph the algorithm runs on. |
661 NetworkSimplex(const GR& graph) : |
691 NetworkSimplex(const GR& graph) : |
662 _graph(graph), |
692 _graph(graph), |
663 _plower(NULL), _pupper(NULL), _pcost(NULL), |
693 _plower(NULL), _pupper(NULL), _pcost(NULL), |
664 _psupply(NULL), _pstsup(false), _ptype(GEQ), |
694 _psupply(NULL), _pstsup(false), _stype(GEQ), |
665 _flow_map(NULL), _potential_map(NULL), |
695 _flow_map(NULL), _potential_map(NULL), |
666 _local_flow(false), _local_potential(false), |
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 && |
702 // Check the value types |
670 std::numeric_limits<Flow>::is_signed, |
703 LEMON_ASSERT(std::numeric_limits<Flow>::is_signed, |
671 "The flow type of NetworkSimplex must be signed integer"); |
704 "The flow type of NetworkSimplex must be signed"); |
672 LEMON_ASSERT(std::numeric_limits<Cost>::is_integer && |
705 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, |
673 std::numeric_limits<Cost>::is_signed, |
706 "The cost type of NetworkSimplex must be signed"); |
674 "The cost type of NetworkSimplex must be signed integer"); |
|
675 } |
707 } |
676 |
708 |
677 /// Destructor. |
709 /// Destructor. |
678 ~NetworkSimplex() { |
710 ~NetworkSimplex() { |
679 if (_local_flow) delete _flow_map; |
711 if (_local_flow) delete _flow_map; |
687 /// @{ |
719 /// @{ |
688 |
720 |
689 /// \brief Set the lower bounds on the arcs. |
721 /// \brief Set the lower bounds on the arcs. |
690 /// |
722 /// |
691 /// This function sets the lower bounds on the arcs. |
723 /// This function sets the lower bounds on the arcs. |
692 /// If neither this function nor \ref boundMaps() is used before |
724 /// If it is not used before calling \ref run(), the lower bounds |
693 /// calling \ref run(), the lower bounds will be set to zero |
725 /// will be set to zero on all arcs. |
694 /// on all arcs. |
|
695 /// |
726 /// |
696 /// \param map An arc map storing the lower bounds. |
727 /// \param map An arc map storing the lower bounds. |
697 /// Its \c Value type must be convertible to the \c Flow type |
728 /// Its \c Value type must be convertible to the \c Flow type |
698 /// of the algorithm. |
729 /// of the algorithm. |
699 /// |
730 /// |
700 /// \return <tt>(*this)</tt> |
731 /// \return <tt>(*this)</tt> |
701 template <typename LOWER> |
732 template <typename LowerMap> |
702 NetworkSimplex& lowerMap(const LOWER& map) { |
733 NetworkSimplex& lowerMap(const LowerMap& map) { |
703 delete _plower; |
734 delete _plower; |
704 _plower = new FlowArcMap(_graph); |
735 _plower = new FlowArcMap(_graph); |
705 for (ArcIt a(_graph); a != INVALID; ++a) { |
736 for (ArcIt a(_graph); a != INVALID; ++a) { |
706 (*_plower)[a] = map[a]; |
737 (*_plower)[a] = map[a]; |
707 } |
738 } |
709 } |
740 } |
710 |
741 |
711 /// \brief Set the upper bounds (capacities) on the arcs. |
742 /// \brief Set the upper bounds (capacities) on the arcs. |
712 /// |
743 /// |
713 /// This function sets the upper bounds (capacities) on the arcs. |
744 /// This function sets the upper bounds (capacities) on the arcs. |
714 /// If none of the functions \ref upperMap(), \ref capacityMap() |
745 /// If it is not used before calling \ref run(), the upper bounds |
715 /// and \ref boundMaps() is used before calling \ref run(), |
746 /// will be set to \ref INF on all arcs (i.e. the flow value will be |
716 /// the upper bounds (capacities) will be set to |
747 /// unbounded from above on each arc). |
717 /// \c std::numeric_limits<Flow>::max() on all arcs. |
|
718 /// |
748 /// |
719 /// \param map An arc map storing the upper bounds. |
749 /// \param map An arc map storing the upper bounds. |
720 /// Its \c Value type must be convertible to the \c Flow type |
750 /// Its \c Value type must be convertible to the \c Flow type |
721 /// of the algorithm. |
751 /// of the algorithm. |
722 /// |
752 /// |
723 /// \return <tt>(*this)</tt> |
753 /// \return <tt>(*this)</tt> |
724 template<typename UPPER> |
754 template<typename UpperMap> |
725 NetworkSimplex& upperMap(const UPPER& map) { |
755 NetworkSimplex& upperMap(const UpperMap& map) { |
726 delete _pupper; |
756 delete _pupper; |
727 _pupper = new FlowArcMap(_graph); |
757 _pupper = new FlowArcMap(_graph); |
728 for (ArcIt a(_graph); a != INVALID; ++a) { |
758 for (ArcIt a(_graph); a != INVALID; ++a) { |
729 (*_pupper)[a] = map[a]; |
759 (*_pupper)[a] = map[a]; |
730 } |
760 } |
731 return *this; |
761 return *this; |
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 before |
|
749 /// calling \ref run(), the lower bounds will be set to zero |
|
750 /// 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 to |
|
754 /// \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 the |
|
760 /// \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 /// \brief Set the costs of the arcs. |
764 /// \brief Set the costs of the arcs. |
772 /// |
765 /// |
773 /// This function sets the costs of the arcs. |
766 /// This function sets the costs of the arcs. |
774 /// If it is not used before calling \ref run(), the costs |
767 /// If it is not used before calling \ref run(), the costs |
775 /// will be set to \c 1 on all arcs. |
768 /// will be set to \c 1 on all arcs. |
817 /// This function sets a single source node and a single target node |
810 /// This function sets a single source node and a single target node |
818 /// and the required flow value. |
811 /// and the required flow value. |
819 /// If neither this function nor \ref supplyMap() is used before |
812 /// If neither this function nor \ref supplyMap() is used before |
820 /// calling \ref run(), the supply of each node will be set to zero. |
813 /// calling \ref run(), the supply of each node will be set to zero. |
821 /// (It makes sense only if non-zero lower bounds are given.) |
814 /// (It makes sense only if non-zero 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 /// \param s The source node. |
820 /// \param s The source node. |
824 /// \param t The target node. |
821 /// \param t The target node. |
825 /// \param k The required amount of flow from node \c s to node \c t |
822 /// \param k The required amount of flow from node \c s to node \c t |
826 /// (i.e. the supply of \c s and the demand of \c t). |
823 /// (i.e. the supply of \c s and the demand of \c t). |
894 |
891 |
895 /// \brief Run the algorithm. |
892 /// \brief Run the algorithm. |
896 /// |
893 /// |
897 /// This function runs the algorithm. |
894 /// This function runs the algorithm. |
898 /// The paramters can be specified using functions \ref lowerMap(), |
895 /// The paramters can be specified using functions \ref lowerMap(), |
899 /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(), |
896 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), |
900 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), |
897 /// \ref supplyType(), \ref flowMap() and \ref potentialMap(). |
901 /// \ref problemType(), \ref flowMap() and \ref potentialMap(). |
|
902 /// For example, |
898 /// For example, |
903 /// \code |
899 /// \code |
904 /// NetworkSimplex<ListDigraph> ns(graph); |
900 /// NetworkSimplex<ListDigraph> ns(graph); |
905 /// ns.boundMaps(lower, upper).costMap(cost) |
901 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
906 /// .supplyMap(sup).run(); |
902 /// .supplyMap(sup).run(); |
907 /// \endcode |
903 /// \endcode |
908 /// |
904 /// |
909 /// This function can be called more than once. All the parameters |
905 /// This function can be called more than once. All the parameters |
910 /// that have been given are kept for the next call, unless |
906 /// that have been given are kept for the next call, unless |
912 /// have to be set again. See \ref reset() for examples. |
908 /// have to be set again. See \ref reset() for examples. |
913 /// |
909 /// |
914 /// \param pivot_rule The pivot rule that will be used during the |
910 /// \param pivot_rule The pivot rule that will be used during the |
915 /// algorithm. For more information see \ref PivotRule. |
911 /// algorithm. For more information see \ref PivotRule. |
916 /// |
912 /// |
917 /// \return \c true if a feasible flow can be found. |
913 /// \return \c INFEASIBLE if no feasible flow exists, |
918 bool run(PivotRule pivot_rule = BLOCK_SEARCH) { |
914 /// \n \c OPTIMAL if the problem has optimal solution |
919 return init() && start(pivot_rule); |
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 |
922 /// \brief Reset all the parameters that have been given before. |
927 /// \brief Reset all the parameters that have been given before. |
923 /// |
928 /// |
924 /// This function resets all the paramaters that have been given |
929 /// This function resets all the paramaters that have been given |
925 /// before using functions \ref lowerMap(), \ref upperMap(), |
930 /// before using functions \ref lowerMap(), \ref upperMap(), |
926 /// \ref capacityMap(), \ref boundMaps(), \ref costMap(), |
931 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(), |
927 /// \ref supplyMap(), \ref stSupply(), \ref problemType(), |
|
928 /// \ref flowMap() and \ref potentialMap(). |
932 /// \ref flowMap() and \ref potentialMap(). |
929 /// |
933 /// |
930 /// It is useful for multiple run() calls. If this function is not |
934 /// It is useful for multiple run() calls. If this function is not |
931 /// used, all the parameters given before are kept for the next |
935 /// used, all the parameters given before are kept for the next |
932 /// \ref run() call. |
936 /// \ref run() call. |
934 /// For example, |
938 /// For example, |
935 /// \code |
939 /// \code |
936 /// NetworkSimplex<ListDigraph> ns(graph); |
940 /// NetworkSimplex<ListDigraph> ns(graph); |
937 /// |
941 /// |
938 /// // First run |
942 /// // First run |
939 /// ns.lowerMap(lower).capacityMap(cap).costMap(cost) |
943 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) |
940 /// .supplyMap(sup).run(); |
944 /// .supplyMap(sup).run(); |
941 /// |
945 /// |
942 /// // Run again with modified cost map (reset() is not called, |
946 /// // Run again with modified cost map (reset() is not called, |
943 /// // so only the cost map have to be set again) |
947 /// // so only the cost map have to be set again) |
944 /// cost[e] += 100; |
948 /// cost[e] += 100; |
945 /// ns.costMap(cost).run(); |
949 /// ns.costMap(cost).run(); |
946 /// |
950 /// |
947 /// // Run again from scratch using reset() |
951 /// // Run again from scratch using reset() |
948 /// // (the lower bounds will be set to zero on all arcs) |
952 /// // (the lower bounds will be set to zero on all arcs) |
949 /// ns.reset(); |
953 /// ns.reset(); |
950 /// ns.capacityMap(cap).costMap(cost) |
954 /// ns.upperMap(capacity).costMap(cost) |
951 /// .supplyMap(sup).run(); |
955 /// .supplyMap(sup).run(); |
952 /// \endcode |
956 /// \endcode |
953 /// |
957 /// |
954 /// \return <tt>(*this)</tt> |
958 /// \return <tt>(*this)</tt> |
955 NetworkSimplex& reset() { |
959 NetworkSimplex& reset() { |
1099 _last_succ.resize(all_node_num); |
1103 _last_succ.resize(all_node_num); |
1100 _state.resize(all_arc_num); |
1104 _state.resize(all_arc_num); |
1101 |
1105 |
1102 // Initialize node related data |
1106 // Initialize node related data |
1103 bool valid_supply = true; |
1107 bool valid_supply = true; |
1104 Flow sum_supply = 0; |
1108 _sum_supply = 0; |
1105 if (!_pstsup && !_psupply) { |
1109 if (!_pstsup && !_psupply) { |
1106 _pstsup = true; |
1110 _pstsup = true; |
1107 _psource = _ptarget = NodeIt(_graph); |
1111 _psource = _ptarget = NodeIt(_graph); |
1108 _pstflow = 0; |
1112 _pstflow = 0; |
1109 } |
1113 } |
1110 if (_psupply) { |
1114 if (_psupply) { |
1111 int i = 0; |
1115 int i = 0; |
1112 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
1116 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
1113 _node_id[n] = i; |
1117 _node_id[n] = i; |
1114 _supply[i] = (*_psupply)[n]; |
1118 _supply[i] = (*_psupply)[n]; |
1115 sum_supply += _supply[i]; |
1119 _sum_supply += _supply[i]; |
1116 } |
1120 } |
1117 valid_supply = (_ptype == GEQ && sum_supply <= 0) || |
1121 valid_supply = (_stype == GEQ && _sum_supply <= 0) || |
1118 (_ptype == LEQ && sum_supply >= 0); |
1122 (_stype == LEQ && _sum_supply >= 0); |
1119 } else { |
1123 } else { |
1120 int i = 0; |
1124 int i = 0; |
1121 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
1125 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { |
1122 _node_id[n] = i; |
1126 _node_id[n] = i; |
1123 _supply[i] = 0; |
1127 _supply[i] = 0; |
1125 _supply[_node_id[_psource]] = _pstflow; |
1129 _supply[_node_id[_psource]] = _pstflow; |
1126 _supply[_node_id[_ptarget]] = -_pstflow; |
1130 _supply[_node_id[_ptarget]] = -_pstflow; |
1127 } |
1131 } |
1128 if (!valid_supply) return false; |
1132 if (!valid_supply) return false; |
1129 |
1133 |
1130 // Infinite capacity value |
|
1131 Flow inf_cap = |
|
1132 std::numeric_limits<Flow>::has_infinity ? |
|
1133 std::numeric_limits<Flow>::infinity() : |
|
1134 std::numeric_limits<Flow>::max(); |
|
1135 |
|
1136 // Initialize artifical cost |
1134 // Initialize artifical cost |
1137 Cost art_cost; |
1135 Cost ART_COST; |
1138 if (std::numeric_limits<Cost>::is_exact) { |
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 } else { |
1138 } else { |
1141 art_cost = std::numeric_limits<Cost>::min(); |
1139 ART_COST = std::numeric_limits<Cost>::min(); |
1142 for (int i = 0; i != _arc_num; ++i) { |
1140 for (int i = 0; i != _arc_num; ++i) { |
1143 if (_cost[i] > art_cost) art_cost = _cost[i]; |
1141 if (_cost[i] > ART_COST) ART_COST = _cost[i]; |
1144 } |
1142 } |
1145 art_cost = (art_cost + 1) * _node_num; |
1143 ART_COST = (ART_COST + 1) * _node_num; |
1146 } |
1144 } |
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; |
|
1215 |
1145 |
1216 // Set data for the artificial root node |
1146 // Set data for the artificial root node |
1217 _root = _node_num; |
1147 _root = _node_num; |
1218 _parent[_root] = -1; |
1148 _parent[_root] = -1; |
1219 _pred[_root] = -1; |
1149 _pred[_root] = -1; |
1220 _thread[_root] = 0; |
1150 _thread[_root] = 0; |
1221 _rev_thread[0] = _root; |
1151 _rev_thread[0] = _root; |
1222 _succ_num[_root] = all_node_num; |
1152 _succ_num[_root] = all_node_num; |
1223 _last_succ[_root] = _root - 1; |
1153 _last_succ[_root] = _root - 1; |
1224 _supply[_root] = -sum_supply; |
1154 _supply[_root] = -_sum_supply; |
1225 if (sum_supply < 0) { |
1155 if (_sum_supply < 0) { |
1226 _pi[_root] = -art_cost; |
1156 _pi[_root] = -ART_COST; |
1227 } else { |
1157 } else { |
1228 _pi[_root] = art_cost; |
1158 _pi[_root] = ART_COST; |
1229 } |
1159 } |
1230 |
1160 |
1231 // Store the arcs in a mixed order |
1161 // Store the arcs in a mixed order |
1232 int k = std::max(int(std::sqrt(double(_arc_num))), 10); |
1162 int k = std::max(int(std::sqrt(double(_arc_num))), 10); |
1233 int i = 0; |
1163 int i = 0; |
1538 case CANDIDATE_LIST: |
1479 case CANDIDATE_LIST: |
1539 return start<CandidateListPivotRule>(); |
1480 return start<CandidateListPivotRule>(); |
1540 case ALTERING_LIST: |
1481 case ALTERING_LIST: |
1541 return start<AlteringListPivotRule>(); |
1482 return start<AlteringListPivotRule>(); |
1542 } |
1483 } |
1543 return false; |
1484 return INFEASIBLE; // avoid warning |
1544 } |
1485 } |
1545 |
1486 |
1546 template <typename PivotRuleImpl> |
1487 template <typename PivotRuleImpl> |
1547 bool start() { |
1488 ProblemType start() { |
1548 PivotRuleImpl pivot(*this); |
1489 PivotRuleImpl pivot(*this); |
1549 |
1490 |
1550 // Execute the Network Simplex algorithm |
1491 // Execute the Network Simplex algorithm |
1551 while (pivot.findEnteringArc()) { |
1492 while (pivot.findEnteringArc()) { |
1552 findJoinNode(); |
1493 findJoinNode(); |
1553 bool change = findLeavingArc(); |
1494 bool change = findLeavingArc(); |
|
1495 if (delta >= INF) return UNBOUNDED; |
1554 changeFlow(change); |
1496 changeFlow(change); |
1555 if (change) { |
1497 if (change) { |
1556 updateTreeStructure(); |
1498 updateTreeStructure(); |
1557 updatePotential(); |
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 } |
1560 |
1519 |
1561 // Copy flow values to _flow_map |
1520 // Copy flow values to _flow_map |
1562 if (_plower) { |
1521 if (_plower) { |