Changes in lemon/network_simplex.h [618:b95898314e09:643:f3792d5bb294] in lemon-1.2
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
lemon/network_simplex.h
r618 r643 31 31 #include <lemon/core.h> 32 32 #include <lemon/math.h> 33 #include <lemon/maps.h>34 #include <lemon/circulation.h>35 #include <lemon/adaptors.h>36 33 37 34 namespace lemon { … … 51 48 /// In general this class is the fastest implementation available 52 49 /// in LEMON for the minimum cost flow problem. 53 /// Moreover it supports both direction of the supply/demand inequality 54 /// constraints. For more information see \ref ProblemType. 50 /// Moreover it supports both directions of the supply/demand inequality 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 58 /// \tparam GR The digraph type the algorithm runs on. 57 /// \tparam FThe value type used for flow amounts, capacity bounds59 /// \tparam V The value type used for flow amounts, capacity bounds 58 60 /// and supply values in the algorithm. By default it is \c int. 59 61 /// \tparam C The value type used for costs and potentials in the 60 /// algorithm. By default it is the same as \c F.62 /// algorithm. By default it is the same as \c V. 61 63 /// 62 64 /// \warning Both value types must be signed and all input data must … … 66 68 /// implementations, from which the most efficient one is used 67 69 /// by default. For more information see \ref PivotRule. 68 template <typename GR, typename F = int, typename C = F>70 template <typename GR, typename V = int, typename C = V> 69 71 class NetworkSimplex 70 72 { 71 73 public: 72 74 73 /// The flow type of the algorithm74 typedef F Flow;75 /// The cost type of the algorithm75 /// The type of the flow amounts, capacity bounds and supply values 76 typedef V Value; 77 /// The type of the arc costs 76 78 typedef C Cost; 77 #ifdef DOXYGEN78 /// The type of the flow map79 typedef GR::ArcMap<Flow> FlowMap;80 /// The type of the potential map81 typedef GR::NodeMap<Cost> PotentialMap;82 #else83 /// The type of the flow map84 typedef typename GR::template ArcMap<Flow> FlowMap;85 /// The type of the potential map86 typedef typename GR::template NodeMap<Cost> PotentialMap;87 #endif88 79 89 80 public: 90 81 91 /// \brief Enum type for selecting the pivot rule. 92 /// 93 /// Enum type for selecting the pivot rule for the \ref run() 94 /// function. 95 /// 96 /// \ref NetworkSimplex provides five different pivot rule 97 /// implementations that significantly affect the running time 98 /// of the algorithm. 99 /// By default \ref BLOCK_SEARCH "Block Search" is used, which 100 /// proved to be the most efficient and the most robust on various 101 /// test inputs according to our benchmark tests. 102 /// However another pivot rule can be selected using the \ref run() 103 /// function with the proper parameter. 104 enum PivotRule { 105 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 82 /// \brief Problem type constants for the \c run() function. 83 /// 84 /// Enum type containing the problem type constants that can be 85 /// returned by the \ref run() function of the algorithm. 86 enum ProblemType { 87 /// The problem has no feasible solution (flow). 88 INFEASIBLE, 89 /// The problem has optimal solution (i.e. it is feasible and 90 /// bounded), and the algorithm has found optimal flow and node 91 /// potentials (primal and dual solutions). 92 OPTIMAL, 93 /// The objective function of the problem is unbounded, i.e. 94 /// there is a directed cycle having negative total cost and 95 /// infinite upper bound. 96 UNBOUNDED 132 97 }; 133 98 134 /// \brief Enum type for selecting the problem type.135 /// 136 /// Enum type for selecting the problem type, i.e. the direction of137 /// the inequalities in the supply/demand constraints of the138 /// \ref min_cost_flow "minimum cost flow problem".139 /// 140 /// The default problemtype is \c GEQ, since this form is supported99 /// \brief Constants for selecting the type of the supply constraints. 100 /// 101 /// Enum type containing constants for selecting the supply type, 102 /// i.e. the direction of the inequalities in the supply/demand 103 /// constraints of the \ref min_cost_flow "minimum cost flow problem". 104 /// 105 /// The default supply type is \c GEQ, since this form is supported 141 106 /// by other minimum cost flow algorithms and the \ref Circulation 142 /// algorithm as well.143 /// The \c LEQ problem type can be selected using the \ref problemType()107 /// algorithm, as well. 108 /// The \c LEQ problem type can be selected using the \ref supplyType() 144 109 /// function. 145 110 /// 146 /// Note that the equality form is a special case of both problem type.147 enum ProblemType {148 149 /// This option means that there are "<em>greater or equal</em>"150 /// constraints in the defintion, i.e. the exact formulation of the151 /// problem is the following.111 /// Note that the equality form is a special case of both supply types. 112 enum SupplyType { 113 114 /// This option means that there are <em>"greater or equal"</em> 115 /// supply/demand constraints in the definition, i.e. the exact 116 /// formulation of the problem is the following. 152 117 /** 153 118 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] … … 165 130 CARRY_SUPPLIES = GEQ, 166 131 167 /// This option means that there are "<em>less or equal</em>"168 /// constraints in the defintion, i.e. the exact formulation of the169 /// problem is the following.132 /// This option means that there are <em>"less or equal"</em> 133 /// supply/demand constraints in the definition, i.e. the exact 134 /// formulation of the problem is the following. 170 135 /** 171 136 \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f] … … 183 148 SATISFY_DEMANDS = LEQ 184 149 }; 185 150 151 /// \brief Constants for selecting the pivot rule. 152 /// 153 /// Enum type containing constants for selecting the pivot rule for 154 /// the \ref run() function. 155 /// 156 /// \ref NetworkSimplex provides five different pivot rule 157 /// implementations that significantly affect the running time 158 /// of the algorithm. 159 /// By default \ref BLOCK_SEARCH "Block Search" is used, which 160 /// proved to be the most efficient and the most robust on various 161 /// test inputs according to our benchmark tests. 162 /// However another pivot rule can be selected using the \ref run() 163 /// function with the proper parameter. 164 enum PivotRule { 165 166 /// The First Eligible pivot rule. 167 /// The next eligible arc is selected in a wraparound fashion 168 /// in every iteration. 169 FIRST_ELIGIBLE, 170 171 /// The Best Eligible pivot rule. 172 /// The best eligible arc is selected in every iteration. 173 BEST_ELIGIBLE, 174 175 /// The Block Search pivot rule. 176 /// A specified number of arcs are examined in every iteration 177 /// in a wraparound fashion and the best eligible arc is selected 178 /// from this block. 179 BLOCK_SEARCH, 180 181 /// The Candidate List pivot rule. 182 /// In a major iteration a candidate list is built from eligible arcs 183 /// in a wraparound fashion and in the following minor iterations 184 /// the best eligible arc is selected from this list. 185 CANDIDATE_LIST, 186 187 /// The Altering Candidate List pivot rule. 188 /// It is a modified version of the Candidate List method. 189 /// It keeps only the several best eligible arcs from the former 190 /// candidate list and extends this list in every iteration. 191 ALTERING_LIST 192 }; 193 186 194 private: 187 195 188 196 TEMPLATE_DIGRAPH_TYPEDEFS(GR); 189 190 typedef typename GR::template ArcMap<Flow> FlowArcMap;191 typedef typename GR::template ArcMap<Cost> CostArcMap;192 typedef typename GR::template NodeMap<Flow> FlowNodeMap;193 197 194 198 typedef std::vector<Arc> ArcVector; … … 196 200 typedef std::vector<int> IntVector; 197 201 typedef std::vector<bool> BoolVector; 198 typedef std::vector< Flow> FlowVector;202 typedef std::vector<Value> ValueVector; 199 203 typedef std::vector<Cost> CostVector; 200 204 … … 214 218 215 219 // Parameters of the problem 216 FlowArcMap *_plower; 217 FlowArcMap *_pupper; 218 CostArcMap *_pcost; 219 FlowNodeMap *_psupply; 220 bool _pstsup; 221 Node _psource, _ptarget; 222 Flow _pstflow; 223 ProblemType _ptype; 224 225 // Result maps 226 FlowMap *_flow_map; 227 PotentialMap *_potential_map; 228 bool _local_flow; 229 bool _local_potential; 220 bool _have_lower; 221 SupplyType _stype; 222 Value _sum_supply; 230 223 231 224 // Data structures for storing the digraph 232 225 IntNodeMap _node_id; 233 ArcVector _arc_ref;226 IntArcMap _arc_id; 234 227 IntVector _source; 235 228 IntVector _target; 236 229 237 230 // Node and arc data 238 FlowVector _cap; 231 ValueVector _lower; 232 ValueVector _upper; 233 ValueVector _cap; 239 234 CostVector _cost; 240 FlowVector _supply;241 FlowVector _flow;235 ValueVector _supply; 236 ValueVector _flow; 242 237 CostVector _pi; 243 238 … … 258 253 int first, second, right, last; 259 254 int stem, par_stem, new_stem; 260 Flow delta; 255 Value delta; 256 257 public: 258 259 /// \brief Constant for infinite upper bounds (capacities). 260 /// 261 /// Constant for infinite upper bounds (capacities). 262 /// It is \c std::numeric_limits<Value>::infinity() if available, 263 /// \c std::numeric_limits<Value>::max() otherwise. 264 const Value INF; 261 265 262 266 private: … … 660 664 /// \param graph The digraph the algorithm runs on. 661 665 NetworkSimplex(const GR& graph) : 662 _graph(graph), 663 _plower(NULL), _pupper(NULL), _pcost(NULL), 664 _psupply(NULL), _pstsup(false), _ptype(GEQ), 665 _flow_map(NULL), _potential_map(NULL), 666 _local_flow(false), _local_potential(false), 667 _node_id(graph) 666 _graph(graph), _node_id(graph), _arc_id(graph), 667 INF(std::numeric_limits<Value>::has_infinity ? 668 std::numeric_limits<Value>::infinity() : 669 std::numeric_limits<Value>::max()) 668 670 { 669 LEMON_ASSERT(std::numeric_limits<Flow>::is_integer && 670 std::numeric_limits<Flow>::is_signed, 671 "The flow type of NetworkSimplex must be signed integer"); 672 LEMON_ASSERT(std::numeric_limits<Cost>::is_integer && 673 std::numeric_limits<Cost>::is_signed, 674 "The cost type of NetworkSimplex must be signed integer"); 675 } 676 677 /// Destructor. 678 ~NetworkSimplex() { 679 if (_local_flow) delete _flow_map; 680 if (_local_potential) delete _potential_map; 671 // Check the value types 672 LEMON_ASSERT(std::numeric_limits<Value>::is_signed, 673 "The flow type of NetworkSimplex must be signed"); 674 LEMON_ASSERT(std::numeric_limits<Cost>::is_signed, 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; 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; 681 728 } 682 729 … … 690 737 /// 691 738 /// This function sets the lower bounds on the arcs. 692 /// If neither this function nor \ref boundMaps() is used before 693 /// calling \ref run(), the lower bounds will be set to zero 694 /// on all arcs. 739 /// If it is not used before calling \ref run(), the lower bounds 740 /// will be set to zero on all arcs. 695 741 /// 696 742 /// \param map An arc map storing the lower bounds. 697 /// Its \c Value type must be convertible to the \c Flowtype743 /// Its \c Value type must be convertible to the \c Value type 698 744 /// of the algorithm. 699 745 /// 700 746 /// \return <tt>(*this)</tt> 701 template <typename LOWER> 702 NetworkSimplex& lowerMap(const LOWER& map) { 703 delete _plower; 704 _plower = new FlowArcMap(_graph); 747 template <typename LowerMap> 748 NetworkSimplex& lowerMap(const LowerMap& map) { 749 _have_lower = true; 705 750 for (ArcIt a(_graph); a != INVALID; ++a) { 706 (*_plower)[a] = map[a];751 _lower[_arc_id[a]] = map[a]; 707 752 } 708 753 return *this; … … 712 757 /// 713 758 /// This function sets the upper bounds (capacities) on the arcs. 714 /// If none of the functions \ref upperMap(), \ref capacityMap() 715 /// and \ref boundMaps() is used before calling \ref run(), 716 /// the upper bounds (capacities) will be set to 717 /// \c std::numeric_limits<Flow>::max() on all arcs. 759 /// If it is not used before calling \ref run(), the upper bounds 760 /// will be set to \ref INF on all arcs (i.e. the flow value will be 761 /// unbounded from above on each arc). 718 762 /// 719 763 /// \param map An arc map storing the upper bounds. 720 /// Its \c Value type must be convertible to the \c Flowtype764 /// Its \c Value type must be convertible to the \c Value type 721 765 /// of the algorithm. 722 766 /// 723 767 /// \return <tt>(*this)</tt> 724 template<typename UPPER> 725 NetworkSimplex& upperMap(const UPPER& map) { 726 delete _pupper; 727 _pupper = new FlowArcMap(_graph); 768 template<typename UpperMap> 769 NetworkSimplex& upperMap(const UpperMap& map) { 728 770 for (ArcIt a(_graph); a != INVALID; ++a) { 729 (*_pupper)[a] = map[a];771 _upper[_arc_id[a]] = map[a]; 730 772 } 731 773 return *this; 732 }733 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 before749 /// calling \ref run(), the lower bounds will be set to zero750 /// 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 to754 /// \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 the760 /// \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 774 } 770 775 … … 780 785 /// 781 786 /// \return <tt>(*this)</tt> 782 template<typename COST> 783 NetworkSimplex& costMap(const COST& map) { 784 delete _pcost; 785 _pcost = new CostArcMap(_graph); 787 template<typename CostMap> 788 NetworkSimplex& costMap(const CostMap& map) { 786 789 for (ArcIt a(_graph); a != INVALID; ++a) { 787 (*_pcost)[a] = map[a];790 _cost[_arc_id[a]] = map[a]; 788 791 } 789 792 return *this; … … 798 801 /// 799 802 /// \param map A node map storing the supply values. 800 /// Its \c Value type must be convertible to the \c Flowtype803 /// Its \c Value type must be convertible to the \c Value type 801 804 /// of the algorithm. 802 805 /// 803 806 /// \return <tt>(*this)</tt> 804 template<typename SUP> 805 NetworkSimplex& supplyMap(const SUP& map) { 806 delete _psupply; 807 _pstsup = false; 808 _psupply = new FlowNodeMap(_graph); 807 template<typename SupplyMap> 808 NetworkSimplex& supplyMap(const SupplyMap& map) { 809 809 for (NodeIt n(_graph); n != INVALID; ++n) { 810 (*_psupply)[n] = map[n];810 _supply[_node_id[n]] = map[n]; 811 811 } 812 812 return *this; … … 821 821 /// (It makes sense only if non-zero lower bounds are given.) 822 822 /// 823 /// Using this function has the same effect as using \ref supplyMap() 824 /// with such a map in which \c k is assigned to \c s, \c -k is 825 /// assigned to \c t and all other nodes have zero supply value. 826 /// 823 827 /// \param s The source node. 824 828 /// \param t The target node. … … 827 831 /// 828 832 /// \return <tt>(*this)</tt> 829 NetworkSimplex& stSupply(const Node& s, const Node& t, Flow k) { 830 delete _psupply; 831 _psupply = NULL; 832 _pstsup = true; 833 _psource = s; 834 _ptarget = t; 835 _pstflow = k; 833 NetworkSimplex& stSupply(const Node& s, const Node& t, Value 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; 836 839 return *this; 837 840 } 838 841 839 /// \brief Set the problem type.840 /// 841 /// This function sets the problem type for the algorithm.842 /// If it is not used before calling \ref run(), the \ref GEQ problem842 /// \brief Set the type of the supply constraints. 843 /// 844 /// This function sets the type of the supply/demand constraints. 845 /// If it is not used before calling \ref run(), the \ref GEQ supply 843 846 /// type will be used. 844 847 /// 845 /// For more information see \ref ProblemType.848 /// For more information see \ref SupplyType. 846 849 /// 847 850 /// \return <tt>(*this)</tt> 848 NetworkSimplex& problemType(ProblemType problem_type) {849 _ ptype = problem_type;851 NetworkSimplex& supplyType(SupplyType supply_type) { 852 _stype = supply_type; 850 853 return *this; 851 854 } 852 855 853 /// \brief Set the flow map.854 ///855 /// This function sets the flow map.856 /// If it is not used before calling \ref run(), an instance will857 /// be allocated automatically. The destructor deallocates this858 /// automatically allocated map, of course.859 ///860 /// \return <tt>(*this)</tt>861 NetworkSimplex& flowMap(FlowMap& map) {862 if (_local_flow) {863 delete _flow_map;864 _local_flow = false;865 }866 _flow_map = ↦867 return *this;868 }869 870 /// \brief Set the potential map.871 ///872 /// This function sets the potential map, which is used for storing873 /// the dual solution.874 /// If it is not used before calling \ref run(), an instance will875 /// be allocated automatically. The destructor deallocates this876 /// automatically allocated map, of course.877 ///878 /// \return <tt>(*this)</tt>879 NetworkSimplex& potentialMap(PotentialMap& map) {880 if (_local_potential) {881 delete _potential_map;882 _local_potential = false;883 }884 _potential_map = ↦885 return *this;886 }887 888 856 /// @} 889 857 … … 897 865 /// This function runs the algorithm. 898 866 /// The paramters can be specified using functions \ref lowerMap(), 899 /// \ref upperMap(), \ref capacityMap(), \ref boundMaps(), 900 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), 901 /// \ref problemType(), \ref flowMap() and \ref potentialMap(). 867 /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 868 /// \ref supplyType(). 902 869 /// For example, 903 870 /// \code 904 871 /// NetworkSimplex<ListDigraph> ns(graph); 905 /// ns. boundMaps(lower,upper).costMap(cost)872 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 906 873 /// .supplyMap(sup).run(); 907 874 /// \endcode … … 911 878 /// \ref reset() is called, thus only the modified parameters 912 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. 913 882 /// 914 883 /// \param pivot_rule The pivot rule that will be used during the 915 884 /// algorithm. For more information see \ref PivotRule. 916 885 /// 917 /// \return \c true if a feasible flow can be found. 918 bool run(PivotRule pivot_rule = BLOCK_SEARCH) { 919 return init() && start(pivot_rule); 886 /// \return \c INFEASIBLE if no feasible flow exists, 887 /// \n \c OPTIMAL if the problem has optimal solution 888 /// (i.e. it is feasible and bounded), and the algorithm has found 889 /// optimal flow and node potentials (primal and dual solutions), 890 /// \n \c UNBOUNDED if the objective function of the problem is 891 /// unbounded, i.e. there is a directed cycle having negative total 892 /// cost and infinite upper bound. 893 /// 894 /// \see ProblemType, PivotRule 895 ProblemType run(PivotRule pivot_rule = BLOCK_SEARCH) { 896 if (!init()) return INFEASIBLE; 897 return start(pivot_rule); 920 898 } 921 899 … … 924 902 /// This function resets all the paramaters that have been given 925 903 /// before using functions \ref lowerMap(), \ref upperMap(), 926 /// \ref capacityMap(), \ref boundMaps(), \ref costMap(), 927 /// \ref supplyMap(), \ref stSupply(), \ref problemType(), 928 /// \ref flowMap() and \ref potentialMap(). 904 /// \ref costMap(), \ref supplyMap(), \ref stSupply(), \ref supplyType(). 929 905 /// 930 906 /// It is useful for multiple run() calls. If this function is not 931 907 /// used, all the parameters given before are kept for the next 932 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. 933 911 /// 934 912 /// For example, … … 937 915 /// 938 916 /// // First run 939 /// ns.lowerMap(lower). capacityMap(cap).costMap(cost)917 /// ns.lowerMap(lower).upperMap(upper).costMap(cost) 940 918 /// .supplyMap(sup).run(); 941 919 /// … … 948 926 /// // (the lower bounds will be set to zero on all arcs) 949 927 /// ns.reset(); 950 /// ns. capacityMap(cap).costMap(cost)928 /// ns.upperMap(capacity).costMap(cost) 951 929 /// .supplyMap(sup).run(); 952 930 /// \endcode … … 954 932 /// \return <tt>(*this)</tt> 955 933 NetworkSimplex& reset() { 956 delete _plower; 957 delete _pupper; 958 delete _pcost; 959 delete _psupply; 960 _plower = NULL; 961 _pupper = NULL; 962 _pcost = NULL; 963 _psupply = NULL; 964 _pstsup = false; 965 _ptype = GEQ; 966 if (_local_flow) delete _flow_map; 967 if (_local_potential) delete _potential_map; 968 _flow_map = NULL; 969 _potential_map = NULL; 970 _local_flow = false; 971 _local_potential = false; 972 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; 943 _stype = GEQ; 973 944 return *this; 974 945 } … … 986 957 /// 987 958 /// This function returns the total cost of the found flow. 988 /// The complexity of the functionis O(e).959 /// Its complexity is O(e). 989 960 /// 990 961 /// \note The return type of the function can be specified as a … … 998 969 /// 999 970 /// \pre \ref run() must be called before using this function. 1000 template <typename Num> 1001 Num totalCost() const { 1002 Num c = 0; 1003 if (_pcost) { 1004 for (ArcIt e(_graph); e != INVALID; ++e) 1005 c += (*_flow_map)[e] * (*_pcost)[e]; 1006 } else { 1007 for (ArcIt e(_graph); e != INVALID; ++e) 1008 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]); 1009 977 } 1010 978 return c; … … 1022 990 /// 1023 991 /// \pre \ref run() must be called before using this function. 1024 Flow flow(const Arc& a) const { 1025 return (*_flow_map)[a]; 1026 } 1027 1028 /// \brief Return a const reference to the flow map. 1029 /// 1030 /// This function returns a const reference to an arc map storing 1031 /// the found flow. 992 Value flow(const Arc& a) const { 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. 1032 1001 /// 1033 1002 /// \pre \ref run() must be called before using this function. 1034 const FlowMap& flowMap() const { 1035 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 } 1036 1008 } 1037 1009 … … 1043 1015 /// \pre \ref run() must be called before using this function. 1044 1016 Cost potential(const Node& n) const { 1045 return (*_potential_map)[n];1046 } 1047 1048 /// \brief Return a const reference to the potential map1049 /// (the dual solution).1050 /// 1051 /// This function returns a const reference to a node map storing1052 /// the found potentials, which form the dual solution ofthe1053 /// \ 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. 1054 1026 /// 1055 1027 /// \pre \ref run() must be called before using this function. 1056 const PotentialMap& potentialMap() const { 1057 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 } 1058 1033 } 1059 1034 … … 1064 1039 // Initialize internal data structures 1065 1040 bool init() { 1066 // Initialize result maps1067 if (!_flow_map) {1068 _flow_map = new FlowMap(_graph);1069 _local_flow = true;1070 }1071 if (!_potential_map) {1072 _potential_map = new PotentialMap(_graph);1073 _local_potential = true;1074 }1075 1076 // Initialize vectors1077 _node_num = countNodes(_graph);1078 _arc_num = countArcs(_graph);1079 int all_node_num = _node_num + 1;1080 int all_arc_num = _arc_num + _node_num;1081 1041 if (_node_num == 0) return false; 1082 1042 1083 _arc_ref.resize(_arc_num); 1084 _source.resize(all_arc_num); 1085 _target.resize(all_arc_num); 1086 1087 _cap.resize(all_arc_num); 1088 _cost.resize(all_arc_num); 1089 _supply.resize(all_node_num); 1090 _flow.resize(all_arc_num); 1091 _pi.resize(all_node_num); 1092 1093 _parent.resize(all_node_num); 1094 _pred.resize(all_node_num); 1095 _forward.resize(all_node_num); 1096 _thread.resize(all_node_num); 1097 _rev_thread.resize(all_node_num); 1098 _succ_num.resize(all_node_num); 1099 _last_succ.resize(all_node_num); 1100 _state.resize(all_arc_num); 1101 1102 // Initialize node related data 1103 bool valid_supply = true; 1104 Flow sum_supply = 0; 1105 if (!_pstsup && !_psupply) { 1106 _pstsup = true; 1107 _psource = _ptarget = NodeIt(_graph); 1108 _pstflow = 0; 1109 } 1110 if (_psupply) { 1111 int i = 0; 1112 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 1113 _node_id[n] = i; 1114 _supply[i] = (*_psupply)[n]; 1115 sum_supply += _supply[i]; 1116 } 1117 valid_supply = (_ptype == GEQ && sum_supply <= 0) || 1118 (_ptype == LEQ && sum_supply >= 0); 1043 // Check the sum of supply values 1044 _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 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; 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 } 1119 1063 } else { 1120 int i = 0; 1121 for (NodeIt n(_graph); n != INVALID; ++n, ++i) { 1122 _node_id[n] = i; 1123 _supply[i] = 0; 1124 } 1125 _supply[_node_id[_psource]] = _pstflow; 1126 _supply[_node_id[_ptarget]] = -_pstflow; 1127 } 1128 if (!valid_supply) return false; 1129 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(); 1064 for (int i = 0; i != _arc_num; ++i) { 1065 _cap[i] = _upper[i]; 1066 } 1067 } 1135 1068 1136 1069 // Initialize artifical cost 1137 Cost art_cost;1070 Cost ART_COST; 1138 1071 if (std::numeric_limits<Cost>::is_exact) { 1139 art_cost= std::numeric_limits<Cost>::max() / 4 + 1;1072 ART_COST = std::numeric_limits<Cost>::max() / 4 + 1; 1140 1073 } else { 1141 art_cost= std::numeric_limits<Cost>::min();1074 ART_COST = std::numeric_limits<Cost>::min(); 1142 1075 for (int i = 0; i != _arc_num; ++i) { 1143 if (_cost[i] > art_cost) art_cost = _cost[i]; 1144 } 1145 art_cost = (art_cost + 1) * _node_num; 1146 } 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 1076 if (_cost[i] > ART_COST) ART_COST = _cost[i]; 1077 } 1078 ART_COST = (ART_COST + 1) * _node_num; 1079 } 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 1216 1087 // Set data for the artificial root node 1217 1088 _root = _node_num; … … 1220 1091 _thread[_root] = 0; 1221 1092 _rev_thread[0] = _root; 1222 _succ_num[_root] = all_node_num;1093 _succ_num[_root] = _node_num + 1; 1223 1094 _last_succ[_root] = _root - 1; 1224 _supply[_root] = -sum_supply; 1225 if (sum_supply < 0) { 1226 _pi[_root] = -art_cost; 1227 } else { 1228 _pi[_root] = art_cost; 1229 } 1230 1231 // Store the arcs in a mixed order 1232 int k = std::max(int(std::sqrt(double(_arc_num))), 10); 1233 int i = 0; 1234 for (ArcIt e(_graph); e != INVALID; ++e) { 1235 _arc_ref[i] = e; 1236 if ((i += k) >= _arc_num) i = (i % k) + 1; 1237 } 1238 1239 // Initialize arc maps 1240 if (_pupper && _pcost) { 1241 for (int i = 0; i != _arc_num; ++i) { 1242 Arc e = _arc_ref[i]; 1243 _source[i] = _node_id[_graph.source(e)]; 1244 _target[i] = _node_id[_graph.target(e)]; 1245 _cap[i] = (*_pupper)[e]; 1246 _cost[i] = (*_pcost)[e]; 1247 _flow[i] = 0; 1248 _state[i] = STATE_LOWER; 1249 } 1250 } else { 1251 for (int i = 0; i != _arc_num; ++i) { 1252 Arc e = _arc_ref[i]; 1253 _source[i] = _node_id[_graph.source(e)]; 1254 _target[i] = _node_id[_graph.target(e)]; 1255 _flow[i] = 0; 1256 _state[i] = STATE_LOWER; 1257 } 1258 if (_pupper) { 1259 for (int i = 0; i != _arc_num; ++i) 1260 _cap[i] = (*_pupper)[_arc_ref[i]]; 1261 } else { 1262 for (int i = 0; i != _arc_num; ++i) 1263 _cap[i] = inf_cap; 1264 } 1265 if (_pcost) { 1266 for (int i = 0; i != _arc_num; ++i) 1267 _cost[i] = (*_pcost)[_arc_ref[i]]; 1268 } else { 1269 for (int i = 0; i != _arc_num; ++i) 1270 _cost[i] = 1; 1271 } 1272 } 1273 1274 // Remove non-zero lower bounds 1275 if (_plower) { 1276 for (int i = 0; i != _arc_num; ++i) { 1277 Flow c = (*_plower)[_arc_ref[i]]; 1278 if (c != 0) { 1279 _cap[i] -= c; 1280 _supply[_source[i]] -= c; 1281 _supply[_target[i]] += c; 1282 } 1283 } 1284 } 1095 _supply[_root] = -_sum_supply; 1096 _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST; 1285 1097 1286 1098 // Add artificial arcs and initialize the spanning tree data structure 1287 1099 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1100 _parent[u] = _root; 1101 _pred[u] = e; 1288 1102 _thread[u] = u + 1; 1289 1103 _rev_thread[u + 1] = u; 1290 1104 _succ_num[u] = 1; 1291 1105 _last_succ[u] = u; 1292 _parent[u] = _root; 1293 _pred[u] = e; 1294 _cost[e] = art_cost; 1295 _cap[e] = inf_cap; 1106 _cost[e] = ART_COST; 1107 _cap[e] = INF; 1296 1108 _state[e] = STATE_TREE; 1297 if (_supply[u] > 0 || (_supply[u] == 0 && sum_supply <= 0)) {1109 if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) { 1298 1110 _flow[e] = _supply[u]; 1299 1111 _forward[u] = true; 1300 _pi[u] = - art_cost+ _pi[_root];1112 _pi[u] = -ART_COST + _pi[_root]; 1301 1113 } else { 1302 1114 _flow[e] = -_supply[u]; 1303 1115 _forward[u] = false; 1304 _pi[u] = art_cost+ _pi[_root];1116 _pi[u] = ART_COST + _pi[_root]; 1305 1117 } 1306 1118 } … … 1337 1149 delta = _cap[in_arc]; 1338 1150 int result = 0; 1339 Flowd;1151 Value d; 1340 1152 int e; 1341 1153 … … 1343 1155 for (int u = first; u != join; u = _parent[u]) { 1344 1156 e = _pred[u]; 1345 d = _forward[u] ? _flow[e] : _cap[e] - _flow[e]; 1157 d = _forward[u] ? 1158 _flow[e] : (_cap[e] == INF ? INF : _cap[e] - _flow[e]); 1346 1159 if (d < delta) { 1347 1160 delta = d; … … 1353 1166 for (int u = second; u != join; u = _parent[u]) { 1354 1167 e = _pred[u]; 1355 d = _forward[u] ? _cap[e] - _flow[e] : _flow[e]; 1168 d = _forward[u] ? 1169 (_cap[e] == INF ? INF : _cap[e] - _flow[e]) : _flow[e]; 1356 1170 if (d <= delta) { 1357 1171 delta = d; … … 1375 1189 // Augment along the cycle 1376 1190 if (delta > 0) { 1377 Flowval = _state[in_arc] * delta;1191 Value val = _state[in_arc] * delta; 1378 1192 _flow[in_arc] += val; 1379 1193 for (int u = _source[in_arc]; u != join; u = _parent[u]) { … … 1527 1341 1528 1342 // Execute the algorithm 1529 boolstart(PivotRule pivot_rule) {1343 ProblemType start(PivotRule pivot_rule) { 1530 1344 // Select the pivot rule implementation 1531 1345 switch (pivot_rule) { … … 1541 1355 return start<AlteringListPivotRule>(); 1542 1356 } 1543 return false;1357 return INFEASIBLE; // avoid warning 1544 1358 } 1545 1359 1546 1360 template <typename PivotRuleImpl> 1547 boolstart() {1361 ProblemType start() { 1548 1362 PivotRuleImpl pivot(*this); 1549 1363 … … 1552 1366 findJoinNode(); 1553 1367 bool change = findLeavingArc(); 1368 if (delta >= INF) return UNBOUNDED; 1554 1369 changeFlow(change); 1555 1370 if (change) { … … 1558 1373 } 1559 1374 } 1560 1561 // Copy flow values to _flow_map 1562 if (_plower) { 1375 1376 // Check feasibility 1377 if (_sum_supply < 0) { 1378 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1379 if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE; 1380 } 1381 } 1382 else if (_sum_supply > 0) { 1383 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1384 if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE; 1385 } 1386 } 1387 else { 1388 for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) { 1389 if (_flow[e] != 0) return INFEASIBLE; 1390 } 1391 } 1392 1393 // Transform the solution and the supply map to the original form 1394 if (_have_lower) { 1563 1395 for (int i = 0; i != _arc_num; ++i) { 1564 Arc e = _arc_ref[i]; 1565 _flow_map->set(e, (*_plower)[e] + _flow[i]); 1566 } 1567 } else { 1568 for (int i = 0; i != _arc_num; ++i) { 1569 _flow_map->set(_arc_ref[i], _flow[i]); 1570 } 1571 } 1572 // Copy potential values to _potential_map 1573 for (NodeIt n(_graph); n != INVALID; ++n) { 1574 _potential_map->set(n, _pi[_node_id[n]]); 1575 } 1576 1577 return true; 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 } 1403 } 1404 1405 return OPTIMAL; 1578 1406 } 1579 1407
Note: See TracChangeset
for help on using the changeset viewer.