COIN-OR::LEMON - Graph Library

Changeset 710:8b0df68370a4 in lemon for lemon


Ignore:
Timestamp:
05/12/09 12:06:40 (15 years ago)
Author:
Peter Kovacs <kpeter@…>
Branch:
default
Phase:
public
Message:

Fix the GEQ/LEQ handling in NetworkSimplex? + improve doc (#291)

  • Fix the optimality conditions for the GEQ/LEQ form.
  • Fix the initialization of the algortihm. It ensures correct solutions and it is much faster for the inequality forms.
  • Fix the pivot rules to search all the arcs that have to be allowed to get in the basis.
  • Better block size for the Block Search pivot rule.
  • Improve documentation of the problem and move it to a separate page.
File:
1 edited

Legend:

Unmodified
Added
Removed
  • lemon/network_simplex.h

    r690 r710  
    2020#define LEMON_NETWORK_SIMPLEX_H
    2121
    22 /// \ingroup min_cost_flow
     22/// \ingroup min_cost_flow_algs
    2323///
    2424/// \file
     
    3434namespace lemon {
    3535
    36   /// \addtogroup min_cost_flow
     36  /// \addtogroup min_cost_flow_algs
    3737  /// @{
    3838
     
    103103    /// constraints of the \ref min_cost_flow "minimum cost flow problem".
    104104    ///
    105     /// The default supply type is \c GEQ, since this form is supported
    106     /// by other minimum cost flow algorithms and the \ref Circulation
    107     /// algorithm, as well.
    108     /// The \c LEQ problem type can be selected using the \ref supplyType()
    109     /// function.
    110     ///
    111     /// Note that the equality form is a special case of both supply types.
     105    /// The default supply type is \c GEQ, the \c LEQ type can be
     106    /// selected using \ref supplyType().
     107    /// The equality form is a special case of both supply types.
    112108    enum SupplyType {
    113 
    114109      /// 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.
    117       /**
    118           \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    119           \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
    120               sup(u) \quad \forall u\in V \f]
    121           \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    122       */
    123       /// It means that the total demand must be greater or equal to the
    124       /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
    125       /// negative) and all the supplies have to be carried out from
    126       /// the supply nodes, but there could be demands that are not
    127       /// satisfied.
     110      /// supply/demand constraints in the definition of the problem.
    128111      GEQ,
    129       /// It is just an alias for the \c GEQ option.
    130       CARRY_SUPPLIES = GEQ,
    131 
    132112      /// 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.
    135       /**
    136           \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
    137           \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
    138               sup(u) \quad \forall u\in V \f]
    139           \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
    140       */
    141       /// It means that the total demand must be less or equal to the
    142       /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
    143       /// positive) and all the demands have to be satisfied, but there
    144       /// could be supplies that are not carried out from the supply
    145       /// nodes.
    146       LEQ,
    147       /// It is just an alias for the \c LEQ option.
    148       SATISFY_DEMANDS = LEQ
     113      /// supply/demand constraints in the definition of the problem.
     114      LEQ
    149115    };
    150116   
     
    216182    int _node_num;
    217183    int _arc_num;
     184    int _all_arc_num;
     185    int _search_arc_num;
    218186
    219187    // Parameters of the problem
     
    278246      const CostVector &_pi;
    279247      int &_in_arc;
    280       int _arc_num;
     248      int _search_arc_num;
    281249
    282250      // Pivot rule data
     
    289257        _source(ns._source), _target(ns._target),
    290258        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
    291         _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
     259        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
     260        _next_arc(0)
    292261      {}
    293262
     
    295264      bool findEnteringArc() {
    296265        Cost c;
    297         for (int e = _next_arc; e < _arc_num; ++e) {
     266        for (int e = _next_arc; e < _search_arc_num; ++e) {
    298267          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    299268          if (c < 0) {
     
    329298      const CostVector &_pi;
    330299      int &_in_arc;
    331       int _arc_num;
     300      int _search_arc_num;
    332301
    333302    public:
     
    337306        _source(ns._source), _target(ns._target),
    338307        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
    339         _in_arc(ns.in_arc), _arc_num(ns._arc_num)
     308        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
    340309      {}
    341310
     
    343312      bool findEnteringArc() {
    344313        Cost c, min = 0;
    345         for (int e = 0; e < _arc_num; ++e) {
     314        for (int e = 0; e < _search_arc_num; ++e) {
    346315          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    347316          if (c < min) {
     
    368337      const CostVector &_pi;
    369338      int &_in_arc;
    370       int _arc_num;
     339      int _search_arc_num;
    371340
    372341      // Pivot rule data
     
    380349        _source(ns._source), _target(ns._target),
    381350        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
    382         _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
     351        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
     352        _next_arc(0)
    383353      {
    384354        // The main parameters of the pivot rule
    385         const double BLOCK_SIZE_FACTOR = 2.0;
     355        const double BLOCK_SIZE_FACTOR = 0.5;
    386356        const int MIN_BLOCK_SIZE = 10;
    387357
    388358        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
    389                                     std::sqrt(double(_arc_num))),
     359                                    std::sqrt(double(_search_arc_num))),
    390360                                MIN_BLOCK_SIZE );
    391361      }
     
    396366        int cnt = _block_size;
    397367        int e, min_arc = _next_arc;
    398         for (e = _next_arc; e < _arc_num; ++e) {
     368        for (e = _next_arc; e < _search_arc_num; ++e) {
    399369          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    400370          if (c < min) {
     
    441411      const CostVector &_pi;
    442412      int &_in_arc;
    443       int _arc_num;
     413      int _search_arc_num;
    444414
    445415      // Pivot rule data
     
    455425        _source(ns._source), _target(ns._target),
    456426        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
    457         _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
     427        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
     428        _next_arc(0)
    458429      {
    459430        // The main parameters of the pivot rule
     
    464435
    465436        _list_length = std::max( int(LIST_LENGTH_FACTOR *
    466                                      std::sqrt(double(_arc_num))),
     437                                     std::sqrt(double(_search_arc_num))),
    467438                                 MIN_LIST_LENGTH );
    468439        _minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
     
    501472        min = 0;
    502473        _curr_length = 0;
    503         for (e = _next_arc; e < _arc_num; ++e) {
     474        for (e = _next_arc; e < _search_arc_num; ++e) {
    504475          c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
    505476          if (c < 0) {
     
    547518      const CostVector &_pi;
    548519      int &_in_arc;
    549       int _arc_num;
     520      int _search_arc_num;
    550521
    551522      // Pivot rule data
     
    575546        _source(ns._source), _target(ns._target),
    576547        _cost(ns._cost), _state(ns._state), _pi(ns._pi),
    577         _in_arc(ns.in_arc), _arc_num(ns._arc_num),
    578         _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
     548        _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
     549        _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
    579550      {
    580551        // The main parameters of the pivot rule
     
    585556
    586557        _block_size = std::max( int(BLOCK_SIZE_FACTOR *
    587                                     std::sqrt(double(_arc_num))),
     558                                    std::sqrt(double(_search_arc_num))),
    588559                                MIN_BLOCK_SIZE );
    589560        _head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
     
    611582        int limit = _head_length;
    612583
    613         for (int e = _next_arc; e < _arc_num; ++e) {
     584        for (int e = _next_arc; e < _search_arc_num; ++e) {
    614585          _cand_cost[e] = _state[e] *
    615586            (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
     
    679650      _arc_num = countArcs(_graph);
    680651      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);
     652      int max_arc_num = _arc_num + 2 * _node_num;
     653
     654      _source.resize(max_arc_num);
     655      _target.resize(max_arc_num);
     656
     657      _lower.resize(_arc_num);
     658      _upper.resize(_arc_num);
     659      _cap.resize(max_arc_num);
     660      _cost.resize(max_arc_num);
    690661      _supply.resize(all_node_num);
    691       _flow.resize(all_arc_num);
     662      _flow.resize(max_arc_num);
    692663      _pi.resize(all_node_num);
    693664
     
    699670      _succ_num.resize(all_node_num);
    700671      _last_succ.resize(all_node_num);
    701       _state.resize(all_arc_num);
     672      _state.resize(max_arc_num);
    702673
    703674      // Copy the graph (store the arcs in a mixed order)
     
    10701041      Cost ART_COST;
    10711042      if (std::numeric_limits<Cost>::is_exact) {
    1072         ART_COST = std::numeric_limits<Cost>::max() / 4 + 1;
     1043        ART_COST = std::numeric_limits<Cost>::max() / 2 + 1;
    10731044      } else {
    10741045        ART_COST = std::numeric_limits<Cost>::min();
     
    10941065      _last_succ[_root] = _root - 1;
    10951066      _supply[_root] = -_sum_supply;
    1096       _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
     1067      _pi[_root] = 0;
    10971068
    10981069      // Add artificial arcs and initialize the spanning tree data structure
    1099       for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
    1100         _parent[u] = _root;
    1101         _pred[u] = e;
    1102         _thread[u] = u + 1;
    1103         _rev_thread[u + 1] = u;
    1104         _succ_num[u] = 1;
    1105         _last_succ[u] = u;
    1106         _cost[e] = ART_COST;
    1107         _cap[e] = INF;
    1108         _state[e] = STATE_TREE;
    1109         if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
    1110           _flow[e] = _supply[u];
    1111           _forward[u] = true;
    1112           _pi[u] = -ART_COST + _pi[_root];
    1113         } else {
    1114           _flow[e] = -_supply[u];
    1115           _forward[u] = false;
    1116           _pi[u] = ART_COST + _pi[_root];
    1117         }
     1070      if (_sum_supply == 0) {
     1071        // EQ supply constraints
     1072        _search_arc_num = _arc_num;
     1073        _all_arc_num = _arc_num + _node_num;
     1074        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
     1075          _parent[u] = _root;
     1076          _pred[u] = e;
     1077          _thread[u] = u + 1;
     1078          _rev_thread[u + 1] = u;
     1079          _succ_num[u] = 1;
     1080          _last_succ[u] = u;
     1081          _cap[e] = INF;
     1082          _state[e] = STATE_TREE;
     1083          if (_supply[u] >= 0) {
     1084            _forward[u] = true;
     1085            _pi[u] = 0;
     1086            _source[e] = u;
     1087            _target[e] = _root;
     1088            _flow[e] = _supply[u];
     1089            _cost[e] = 0;
     1090          } else {
     1091            _forward[u] = false;
     1092            _pi[u] = ART_COST;
     1093            _source[e] = _root;
     1094            _target[e] = u;
     1095            _flow[e] = -_supply[u];
     1096            _cost[e] = ART_COST;
     1097          }
     1098        }
     1099      }
     1100      else if (_sum_supply > 0) {
     1101        // LEQ supply constraints
     1102        _search_arc_num = _arc_num + _node_num;
     1103        int f = _arc_num + _node_num;
     1104        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
     1105          _parent[u] = _root;
     1106          _thread[u] = u + 1;
     1107          _rev_thread[u + 1] = u;
     1108          _succ_num[u] = 1;
     1109          _last_succ[u] = u;
     1110          if (_supply[u] >= 0) {
     1111            _forward[u] = true;
     1112            _pi[u] = 0;
     1113            _pred[u] = e;
     1114            _source[e] = u;
     1115            _target[e] = _root;
     1116            _cap[e] = INF;
     1117            _flow[e] = _supply[u];
     1118            _cost[e] = 0;
     1119            _state[e] = STATE_TREE;
     1120          } else {
     1121            _forward[u] = false;
     1122            _pi[u] = ART_COST;
     1123            _pred[u] = f;
     1124            _source[f] = _root;
     1125            _target[f] = u;
     1126            _cap[f] = INF;
     1127            _flow[f] = -_supply[u];
     1128            _cost[f] = ART_COST;
     1129            _state[f] = STATE_TREE;
     1130            _source[e] = u;
     1131            _target[e] = _root;
     1132            _cap[e] = INF;
     1133            _flow[e] = 0;
     1134            _cost[e] = 0;
     1135            _state[e] = STATE_LOWER;
     1136            ++f;
     1137          }
     1138        }
     1139        _all_arc_num = f;
     1140      }
     1141      else {
     1142        // GEQ supply constraints
     1143        _search_arc_num = _arc_num + _node_num;
     1144        int f = _arc_num + _node_num;
     1145        for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
     1146          _parent[u] = _root;
     1147          _thread[u] = u + 1;
     1148          _rev_thread[u + 1] = u;
     1149          _succ_num[u] = 1;
     1150          _last_succ[u] = u;
     1151          if (_supply[u] <= 0) {
     1152            _forward[u] = false;
     1153            _pi[u] = 0;
     1154            _pred[u] = e;
     1155            _source[e] = _root;
     1156            _target[e] = u;
     1157            _cap[e] = INF;
     1158            _flow[e] = -_supply[u];
     1159            _cost[e] = 0;
     1160            _state[e] = STATE_TREE;
     1161          } else {
     1162            _forward[u] = true;
     1163            _pi[u] = -ART_COST;
     1164            _pred[u] = f;
     1165            _source[f] = u;
     1166            _target[f] = _root;
     1167            _cap[f] = INF;
     1168            _flow[f] = _supply[u];
     1169            _state[f] = STATE_TREE;
     1170            _cost[f] = ART_COST;
     1171            _source[e] = _root;
     1172            _target[e] = u;
     1173            _cap[e] = INF;
     1174            _flow[e] = 0;
     1175            _cost[e] = 0;
     1176            _state[e] = STATE_LOWER;
     1177            ++f;
     1178          }
     1179        }
     1180        _all_arc_num = f;
    11181181      }
    11191182
     
    13751438     
    13761439      // 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         }
     1440      for (int e = _search_arc_num; e != _all_arc_num; ++e) {
     1441        if (_flow[e] != 0) return INFEASIBLE;
    13911442      }
    13921443
     
    14021453        }
    14031454      }
     1455     
     1456      // Shift potentials to meet the requirements of the GEQ/LEQ type
     1457      // optimality conditions
     1458      if (_sum_supply == 0) {
     1459        if (_stype == GEQ) {
     1460          Cost max_pot = std::numeric_limits<Cost>::min();
     1461          for (int i = 0; i != _node_num; ++i) {
     1462            if (_pi[i] > max_pot) max_pot = _pi[i];
     1463          }
     1464          if (max_pot > 0) {
     1465            for (int i = 0; i != _node_num; ++i)
     1466              _pi[i] -= max_pot;
     1467          }
     1468        } else {
     1469          Cost min_pot = std::numeric_limits<Cost>::max();
     1470          for (int i = 0; i != _node_num; ++i) {
     1471            if (_pi[i] < min_pot) min_pot = _pi[i];
     1472          }
     1473          if (min_pot < 0) {
     1474            for (int i = 0; i != _node_num; ++i)
     1475              _pi[i] -= min_pot;
     1476          }
     1477        }
     1478      }
    14041479
    14051480      return OPTIMAL;
Note: See TracChangeset for help on using the changeset viewer.