lemon/network_simplex.h
changeset 880 38213abd2911
parent 862 b6f76c95992e
child 892 cbf32bf95954
equal deleted inserted replaced
28:e2acf237e0e0 29:7d65283874b3
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     1 /* -*- mode: C++; indent-tabs-mode: nil; -*-
     2  *
     2  *
     3  * This file is a part of LEMON, a generic C++ optimization library.
     3  * This file is a part of LEMON, a generic C++ optimization library.
     4  *
     4  *
     5  * Copyright (C) 2003-2009
     5  * Copyright (C) 2003-2010
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     6  * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     7  * (Egervary Research Group on Combinatorial Optimization, EGRES).
     8  *
     8  *
     9  * Permission to use, modify and distribute this software is granted
     9  * Permission to use, modify and distribute this software is granted
    10  * provided that this copyright notice appears in all copies. For
    10  * provided that this copyright notice appears in all copies. For
    95       /// The objective function of the problem is unbounded, i.e.
    95       /// The objective function of the problem is unbounded, i.e.
    96       /// there is a directed cycle having negative total cost and
    96       /// there is a directed cycle having negative total cost and
    97       /// infinite upper bound.
    97       /// infinite upper bound.
    98       UNBOUNDED
    98       UNBOUNDED
    99     };
    99     };
   100     
   100 
   101     /// \brief Constants for selecting the type of the supply constraints.
   101     /// \brief Constants for selecting the type of the supply constraints.
   102     ///
   102     ///
   103     /// Enum type containing constants for selecting the supply type,
   103     /// Enum type containing constants for selecting the supply type,
   104     /// i.e. the direction of the inequalities in the supply/demand
   104     /// i.e. the direction of the inequalities in the supply/demand
   105     /// constraints of the \ref min_cost_flow "minimum cost flow problem".
   105     /// constraints of the \ref min_cost_flow "minimum cost flow problem".
   113       GEQ,
   113       GEQ,
   114       /// This option means that there are <em>"less or equal"</em>
   114       /// This option means that there are <em>"less or equal"</em>
   115       /// supply/demand constraints in the definition of the problem.
   115       /// supply/demand constraints in the definition of the problem.
   116       LEQ
   116       LEQ
   117     };
   117     };
   118     
   118 
   119     /// \brief Constants for selecting the pivot rule.
   119     /// \brief Constants for selecting the pivot rule.
   120     ///
   120     ///
   121     /// Enum type containing constants for selecting the pivot rule for
   121     /// Enum type containing constants for selecting the pivot rule for
   122     /// the \ref run() function.
   122     /// the \ref run() function.
   123     ///
   123     ///
   156       /// It is a modified version of the Candidate List method.
   156       /// It is a modified version of the Candidate List method.
   157       /// It keeps only the several best eligible arcs from the former
   157       /// It keeps only the several best eligible arcs from the former
   158       /// candidate list and extends this list in every iteration.
   158       /// candidate list and extends this list in every iteration.
   159       ALTERING_LIST
   159       ALTERING_LIST
   160     };
   160     };
   161     
   161 
   162   private:
   162   private:
   163 
   163 
   164     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   164     TEMPLATE_DIGRAPH_TYPEDEFS(GR);
   165 
   165 
   166     typedef std::vector<int> IntVector;
   166     typedef std::vector<int> IntVector;
   225     // Temporary data used in the current pivot iteration
   225     // Temporary data used in the current pivot iteration
   226     int in_arc, join, u_in, v_in, u_out, v_out;
   226     int in_arc, join, u_in, v_in, u_out, v_out;
   227     int first, second, right, last;
   227     int first, second, right, last;
   228     int stem, par_stem, new_stem;
   228     int stem, par_stem, new_stem;
   229     Value delta;
   229     Value delta;
   230     
   230 
   231     const Value MAX;
   231     const Value MAX;
   232 
   232 
   233   public:
   233   public:
   234   
   234 
   235     /// \brief Constant for infinite upper bounds (capacities).
   235     /// \brief Constant for infinite upper bounds (capacities).
   236     ///
   236     ///
   237     /// Constant for infinite upper bounds (capacities).
   237     /// Constant for infinite upper bounds (capacities).
   238     /// It is \c std::numeric_limits<Value>::infinity() if available,
   238     /// It is \c std::numeric_limits<Value>::infinity() if available,
   239     /// \c std::numeric_limits<Value>::max() otherwise.
   239     /// \c std::numeric_limits<Value>::max() otherwise.
   496             }
   496             }
   497             if (_curr_length == _list_length) goto search_end;
   497             if (_curr_length == _list_length) goto search_end;
   498           }
   498           }
   499         }
   499         }
   500         if (_curr_length == 0) return false;
   500         if (_curr_length == 0) return false;
   501       
   501 
   502       search_end:        
   502       search_end:
   503         _minor_count = 1;
   503         _minor_count = 1;
   504         _next_arc = e;
   504         _next_arc = e;
   505         return true;
   505         return true;
   506       }
   506       }
   507 
   507 
   606             limit = 0;
   606             limit = 0;
   607             cnt = _block_size;
   607             cnt = _block_size;
   608           }
   608           }
   609         }
   609         }
   610         if (_curr_length == 0) return false;
   610         if (_curr_length == 0) return false;
   611         
   611 
   612       search_end:
   612       search_end:
   613 
   613 
   614         // Make heap of the candidate list (approximating a partial sort)
   614         // Make heap of the candidate list (approximating a partial sort)
   615         make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
   615         make_heap( _candidates.begin(), _candidates.begin() + _curr_length,
   616                    _sort_func );
   616                    _sort_func );
   632     ///
   632     ///
   633     /// The constructor of the class.
   633     /// The constructor of the class.
   634     ///
   634     ///
   635     /// \param graph The digraph the algorithm runs on.
   635     /// \param graph The digraph the algorithm runs on.
   636     /// \param arc_mixing Indicate if the arcs have to be stored in a
   636     /// \param arc_mixing Indicate if the arcs have to be stored in a
   637     /// mixed order in the internal data structure. 
   637     /// mixed order in the internal data structure.
   638     /// In special cases, it could lead to better overall performance,
   638     /// In special cases, it could lead to better overall performance,
   639     /// but it is usually slower. Therefore it is disabled by default.
   639     /// but it is usually slower. Therefore it is disabled by default.
   640     NetworkSimplex(const GR& graph, bool arc_mixing = false) :
   640     NetworkSimplex(const GR& graph, bool arc_mixing = false) :
   641       _graph(graph), _node_id(graph), _arc_id(graph),
   641       _graph(graph), _node_id(graph), _arc_id(graph),
   642       _arc_mixing(arc_mixing),
   642       _arc_mixing(arc_mixing),
   647       // Check the number types
   647       // Check the number types
   648       LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
   648       LEMON_ASSERT(std::numeric_limits<Value>::is_signed,
   649         "The flow type of NetworkSimplex must be signed");
   649         "The flow type of NetworkSimplex must be signed");
   650       LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
   650       LEMON_ASSERT(std::numeric_limits<Cost>::is_signed,
   651         "The cost type of NetworkSimplex must be signed");
   651         "The cost type of NetworkSimplex must be signed");
   652         
   652 
   653       // Reset data structures
   653       // Reset data structures
   654       reset();
   654       reset();
   655     }
   655     }
   656 
   656 
   657     /// \name Parameters
   657     /// \name Parameters
   761       }
   761       }
   762       _supply[_node_id[s]] =  k;
   762       _supply[_node_id[s]] =  k;
   763       _supply[_node_id[t]] = -k;
   763       _supply[_node_id[t]] = -k;
   764       return *this;
   764       return *this;
   765     }
   765     }
   766     
   766 
   767     /// \brief Set the type of the supply constraints.
   767     /// \brief Set the type of the supply constraints.
   768     ///
   768     ///
   769     /// This function sets the type of the supply/demand constraints.
   769     /// This function sets the type of the supply/demand constraints.
   770     /// If it is not used before calling \ref run(), the \ref GEQ supply
   770     /// If it is not used before calling \ref run(), the \ref GEQ supply
   771     /// type will be used.
   771     /// type will be used.
   787 
   787 
   788     /// \brief Run the algorithm.
   788     /// \brief Run the algorithm.
   789     ///
   789     ///
   790     /// This function runs the algorithm.
   790     /// This function runs the algorithm.
   791     /// The paramters can be specified using functions \ref lowerMap(),
   791     /// The paramters can be specified using functions \ref lowerMap(),
   792     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(), 
   792     /// \ref upperMap(), \ref costMap(), \ref supplyMap(), \ref stSupply(),
   793     /// \ref supplyType().
   793     /// \ref supplyType().
   794     /// For example,
   794     /// For example,
   795     /// \code
   795     /// \code
   796     ///   NetworkSimplex<ListDigraph> ns(graph);
   796     ///   NetworkSimplex<ListDigraph> ns(graph);
   797     ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
   797     ///   ns.lowerMap(lower).upperMap(upper).costMap(cost)
   942           _arc_id[a] = i;
   942           _arc_id[a] = i;
   943           _source[i] = _node_id[_graph.source(a)];
   943           _source[i] = _node_id[_graph.source(a)];
   944           _target[i] = _node_id[_graph.target(a)];
   944           _target[i] = _node_id[_graph.target(a)];
   945         }
   945         }
   946       }
   946       }
   947       
   947 
   948       // Reset parameters
   948       // Reset parameters
   949       resetParams();
   949       resetParams();
   950       return *this;
   950       return *this;
   951     }
   951     }
   952     
   952 
   953     /// @}
   953     /// @}
   954 
   954 
   955     /// \name Query Functions
   955     /// \name Query Functions
   956     /// The results of the algorithm can be obtained using these
   956     /// The results of the algorithm can be obtained using these
   957     /// functions.\n
   957     /// functions.\n
  1087       // Initialize arc maps
  1087       // Initialize arc maps
  1088       for (int i = 0; i != _arc_num; ++i) {
  1088       for (int i = 0; i != _arc_num; ++i) {
  1089         _flow[i] = 0;
  1089         _flow[i] = 0;
  1090         _state[i] = STATE_LOWER;
  1090         _state[i] = STATE_LOWER;
  1091       }
  1091       }
  1092       
  1092 
  1093       // Set data for the artificial root node
  1093       // Set data for the artificial root node
  1094       _root = _node_num;
  1094       _root = _node_num;
  1095       _parent[_root] = -1;
  1095       _parent[_root] = -1;
  1096       _pred[_root] = -1;
  1096       _pred[_root] = -1;
  1097       _thread[_root] = 0;
  1097       _thread[_root] = 0;
  1261         }
  1261         }
  1262       }
  1262       }
  1263       // Search the cycle along the path form the second node to the root
  1263       // Search the cycle along the path form the second node to the root
  1264       for (int u = second; u != join; u = _parent[u]) {
  1264       for (int u = second; u != join; u = _parent[u]) {
  1265         e = _pred[u];
  1265         e = _pred[u];
  1266         d = _forward[u] ? 
  1266         d = _forward[u] ?
  1267           (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
  1267           (_cap[e] >= MAX ? INF : _cap[e] - _flow[e]) : _flow[e];
  1268         if (d <= delta) {
  1268         if (d <= delta) {
  1269           delta = d;
  1269           delta = d;
  1270           u_out = u;
  1270           u_out = u;
  1271           result = 2;
  1271           result = 2;
  1565         if (change) {
  1565         if (change) {
  1566           updateTreeStructure();
  1566           updateTreeStructure();
  1567           updatePotential();
  1567           updatePotential();
  1568         }
  1568         }
  1569       }
  1569       }
  1570       
  1570 
  1571       // Check feasibility
  1571       // Check feasibility
  1572       for (int e = _search_arc_num; e != _all_arc_num; ++e) {
  1572       for (int e = _search_arc_num; e != _all_arc_num; ++e) {
  1573         if (_flow[e] != 0) return INFEASIBLE;
  1573         if (_flow[e] != 0) return INFEASIBLE;
  1574       }
  1574       }
  1575 
  1575 
  1582             _supply[_source[i]] += c;
  1582             _supply[_source[i]] += c;
  1583             _supply[_target[i]] -= c;
  1583             _supply[_target[i]] -= c;
  1584           }
  1584           }
  1585         }
  1585         }
  1586       }
  1586       }
  1587       
  1587 
  1588       // Shift potentials to meet the requirements of the GEQ/LEQ type
  1588       // Shift potentials to meet the requirements of the GEQ/LEQ type
  1589       // optimality conditions
  1589       // optimality conditions
  1590       if (_sum_supply == 0) {
  1590       if (_sum_supply == 0) {
  1591         if (_stype == GEQ) {
  1591         if (_stype == GEQ) {
  1592           Cost max_pot = std::numeric_limits<Cost>::min();
  1592           Cost max_pot = std::numeric_limits<Cost>::min();