Index: lemon/network_simplex.h
===================================================================
--- lemon/network_simplex.h (revision 690)
+++ lemon/network_simplex.h (revision 710)
@@ -20,5 +20,5 @@
#define LEMON_NETWORK_SIMPLEX_H
-/// \ingroup min_cost_flow
+/// \ingroup min_cost_flow_algs
///
/// \file
@@ -34,5 +34,5 @@
namespace lemon {
- /// \addtogroup min_cost_flow
+ /// \addtogroup min_cost_flow_algs
/// @{
@@ -103,48 +103,14 @@
/// constraints of the \ref min_cost_flow "minimum cost flow problem".
///
- /// The default supply type is \c GEQ, since this form is supported
- /// by other minimum cost flow algorithms and the \ref Circulation
- /// algorithm, as well.
- /// The \c LEQ problem type can be selected using the \ref supplyType()
- /// function.
- ///
- /// Note that the equality form is a special case of both supply types.
+ /// The default supply type is \c GEQ, the \c LEQ type can be
+ /// selected using \ref supplyType().
+ /// The equality form is a special case of both supply types.
enum SupplyType {
-
/// This option means that there are "greater or equal"
- /// supply/demand constraints in the definition, i.e. the exact
- /// formulation of the problem is the following.
- /**
- \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
- \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \geq
- sup(u) \quad \forall u\in V \f]
- \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
- */
- /// It means that the total demand must be greater or equal to the
- /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
- /// negative) and all the supplies have to be carried out from
- /// the supply nodes, but there could be demands that are not
- /// satisfied.
+ /// supply/demand constraints in the definition of the problem.
GEQ,
- /// It is just an alias for the \c GEQ option.
- CARRY_SUPPLIES = GEQ,
-
/// This option means that there are "less or equal"
- /// supply/demand constraints in the definition, i.e. the exact
- /// formulation of the problem is the following.
- /**
- \f[ \min\sum_{uv\in A} f(uv) \cdot cost(uv) \f]
- \f[ \sum_{uv\in A} f(uv) - \sum_{vu\in A} f(vu) \leq
- sup(u) \quad \forall u\in V \f]
- \f[ lower(uv) \leq f(uv) \leq upper(uv) \quad \forall uv\in A \f]
- */
- /// It means that the total demand must be less or equal to the
- /// total supply (i.e. \f$\sum_{u\in V} sup(u)\f$ must be zero or
- /// positive) and all the demands have to be satisfied, but there
- /// could be supplies that are not carried out from the supply
- /// nodes.
- LEQ,
- /// It is just an alias for the \c LEQ option.
- SATISFY_DEMANDS = LEQ
+ /// supply/demand constraints in the definition of the problem.
+ LEQ
};
@@ -216,4 +182,6 @@
int _node_num;
int _arc_num;
+ int _all_arc_num;
+ int _search_arc_num;
// Parameters of the problem
@@ -278,5 +246,5 @@
const CostVector &_pi;
int &_in_arc;
- int _arc_num;
+ int _search_arc_num;
// Pivot rule data
@@ -289,5 +257,6 @@
_source(ns._source), _target(ns._target),
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
- _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
+ _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
+ _next_arc(0)
{}
@@ -295,5 +264,5 @@
bool findEnteringArc() {
Cost c;
- for (int e = _next_arc; e < _arc_num; ++e) {
+ for (int e = _next_arc; e < _search_arc_num; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < 0) {
@@ -329,5 +298,5 @@
const CostVector &_pi;
int &_in_arc;
- int _arc_num;
+ int _search_arc_num;
public:
@@ -337,5 +306,5 @@
_source(ns._source), _target(ns._target),
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
- _in_arc(ns.in_arc), _arc_num(ns._arc_num)
+ _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num)
{}
@@ -343,5 +312,5 @@
bool findEnteringArc() {
Cost c, min = 0;
- for (int e = 0; e < _arc_num; ++e) {
+ for (int e = 0; e < _search_arc_num; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < min) {
@@ -368,5 +337,5 @@
const CostVector &_pi;
int &_in_arc;
- int _arc_num;
+ int _search_arc_num;
// Pivot rule data
@@ -380,12 +349,13 @@
_source(ns._source), _target(ns._target),
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
- _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
+ _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
+ _next_arc(0)
{
// The main parameters of the pivot rule
- const double BLOCK_SIZE_FACTOR = 2.0;
+ const double BLOCK_SIZE_FACTOR = 0.5;
const int MIN_BLOCK_SIZE = 10;
_block_size = std::max( int(BLOCK_SIZE_FACTOR *
- std::sqrt(double(_arc_num))),
+ std::sqrt(double(_search_arc_num))),
MIN_BLOCK_SIZE );
}
@@ -396,5 +366,5 @@
int cnt = _block_size;
int e, min_arc = _next_arc;
- for (e = _next_arc; e < _arc_num; ++e) {
+ for (e = _next_arc; e < _search_arc_num; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < min) {
@@ -441,5 +411,5 @@
const CostVector &_pi;
int &_in_arc;
- int _arc_num;
+ int _search_arc_num;
// Pivot rule data
@@ -455,5 +425,6 @@
_source(ns._source), _target(ns._target),
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
- _in_arc(ns.in_arc), _arc_num(ns._arc_num), _next_arc(0)
+ _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
+ _next_arc(0)
{
// The main parameters of the pivot rule
@@ -464,5 +435,5 @@
_list_length = std::max( int(LIST_LENGTH_FACTOR *
- std::sqrt(double(_arc_num))),
+ std::sqrt(double(_search_arc_num))),
MIN_LIST_LENGTH );
_minor_limit = std::max( int(MINOR_LIMIT_FACTOR * _list_length),
@@ -501,5 +472,5 @@
min = 0;
_curr_length = 0;
- for (e = _next_arc; e < _arc_num; ++e) {
+ for (e = _next_arc; e < _search_arc_num; ++e) {
c = _state[e] * (_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
if (c < 0) {
@@ -547,5 +518,5 @@
const CostVector &_pi;
int &_in_arc;
- int _arc_num;
+ int _search_arc_num;
// Pivot rule data
@@ -575,6 +546,6 @@
_source(ns._source), _target(ns._target),
_cost(ns._cost), _state(ns._state), _pi(ns._pi),
- _in_arc(ns.in_arc), _arc_num(ns._arc_num),
- _next_arc(0), _cand_cost(ns._arc_num), _sort_func(_cand_cost)
+ _in_arc(ns.in_arc), _search_arc_num(ns._search_arc_num),
+ _next_arc(0), _cand_cost(ns._search_arc_num), _sort_func(_cand_cost)
{
// The main parameters of the pivot rule
@@ -585,5 +556,5 @@
_block_size = std::max( int(BLOCK_SIZE_FACTOR *
- std::sqrt(double(_arc_num))),
+ std::sqrt(double(_search_arc_num))),
MIN_BLOCK_SIZE );
_head_length = std::max( int(HEAD_LENGTH_FACTOR * _block_size),
@@ -611,5 +582,5 @@
int limit = _head_length;
- for (int e = _next_arc; e < _arc_num; ++e) {
+ for (int e = _next_arc; e < _search_arc_num; ++e) {
_cand_cost[e] = _state[e] *
(_cost[e] + _pi[_source[e]] - _pi[_target[e]]);
@@ -679,15 +650,15 @@
_arc_num = countArcs(_graph);
int all_node_num = _node_num + 1;
- int all_arc_num = _arc_num + _node_num;
-
- _source.resize(all_arc_num);
- _target.resize(all_arc_num);
-
- _lower.resize(all_arc_num);
- _upper.resize(all_arc_num);
- _cap.resize(all_arc_num);
- _cost.resize(all_arc_num);
+ int max_arc_num = _arc_num + 2 * _node_num;
+
+ _source.resize(max_arc_num);
+ _target.resize(max_arc_num);
+
+ _lower.resize(_arc_num);
+ _upper.resize(_arc_num);
+ _cap.resize(max_arc_num);
+ _cost.resize(max_arc_num);
_supply.resize(all_node_num);
- _flow.resize(all_arc_num);
+ _flow.resize(max_arc_num);
_pi.resize(all_node_num);
@@ -699,5 +670,5 @@
_succ_num.resize(all_node_num);
_last_succ.resize(all_node_num);
- _state.resize(all_arc_num);
+ _state.resize(max_arc_num);
// Copy the graph (store the arcs in a mixed order)
@@ -1070,5 +1041,5 @@
Cost ART_COST;
if (std::numeric_limits::is_exact) {
- ART_COST = std::numeric_limits::max() / 4 + 1;
+ ART_COST = std::numeric_limits::max() / 2 + 1;
} else {
ART_COST = std::numeric_limits::min();
@@ -1094,26 +1065,118 @@
_last_succ[_root] = _root - 1;
_supply[_root] = -_sum_supply;
- _pi[_root] = _sum_supply < 0 ? -ART_COST : ART_COST;
+ _pi[_root] = 0;
// Add artificial arcs and initialize the spanning tree data structure
- for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
- _parent[u] = _root;
- _pred[u] = e;
- _thread[u] = u + 1;
- _rev_thread[u + 1] = u;
- _succ_num[u] = 1;
- _last_succ[u] = u;
- _cost[e] = ART_COST;
- _cap[e] = INF;
- _state[e] = STATE_TREE;
- if (_supply[u] > 0 || (_supply[u] == 0 && _sum_supply <= 0)) {
- _flow[e] = _supply[u];
- _forward[u] = true;
- _pi[u] = -ART_COST + _pi[_root];
- } else {
- _flow[e] = -_supply[u];
- _forward[u] = false;
- _pi[u] = ART_COST + _pi[_root];
- }
+ if (_sum_supply == 0) {
+ // EQ supply constraints
+ _search_arc_num = _arc_num;
+ _all_arc_num = _arc_num + _node_num;
+ for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+ _parent[u] = _root;
+ _pred[u] = e;
+ _thread[u] = u + 1;
+ _rev_thread[u + 1] = u;
+ _succ_num[u] = 1;
+ _last_succ[u] = u;
+ _cap[e] = INF;
+ _state[e] = STATE_TREE;
+ if (_supply[u] >= 0) {
+ _forward[u] = true;
+ _pi[u] = 0;
+ _source[e] = u;
+ _target[e] = _root;
+ _flow[e] = _supply[u];
+ _cost[e] = 0;
+ } else {
+ _forward[u] = false;
+ _pi[u] = ART_COST;
+ _source[e] = _root;
+ _target[e] = u;
+ _flow[e] = -_supply[u];
+ _cost[e] = ART_COST;
+ }
+ }
+ }
+ else if (_sum_supply > 0) {
+ // LEQ supply constraints
+ _search_arc_num = _arc_num + _node_num;
+ int f = _arc_num + _node_num;
+ for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+ _parent[u] = _root;
+ _thread[u] = u + 1;
+ _rev_thread[u + 1] = u;
+ _succ_num[u] = 1;
+ _last_succ[u] = u;
+ if (_supply[u] >= 0) {
+ _forward[u] = true;
+ _pi[u] = 0;
+ _pred[u] = e;
+ _source[e] = u;
+ _target[e] = _root;
+ _cap[e] = INF;
+ _flow[e] = _supply[u];
+ _cost[e] = 0;
+ _state[e] = STATE_TREE;
+ } else {
+ _forward[u] = false;
+ _pi[u] = ART_COST;
+ _pred[u] = f;
+ _source[f] = _root;
+ _target[f] = u;
+ _cap[f] = INF;
+ _flow[f] = -_supply[u];
+ _cost[f] = ART_COST;
+ _state[f] = STATE_TREE;
+ _source[e] = u;
+ _target[e] = _root;
+ _cap[e] = INF;
+ _flow[e] = 0;
+ _cost[e] = 0;
+ _state[e] = STATE_LOWER;
+ ++f;
+ }
+ }
+ _all_arc_num = f;
+ }
+ else {
+ // GEQ supply constraints
+ _search_arc_num = _arc_num + _node_num;
+ int f = _arc_num + _node_num;
+ for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
+ _parent[u] = _root;
+ _thread[u] = u + 1;
+ _rev_thread[u + 1] = u;
+ _succ_num[u] = 1;
+ _last_succ[u] = u;
+ if (_supply[u] <= 0) {
+ _forward[u] = false;
+ _pi[u] = 0;
+ _pred[u] = e;
+ _source[e] = _root;
+ _target[e] = u;
+ _cap[e] = INF;
+ _flow[e] = -_supply[u];
+ _cost[e] = 0;
+ _state[e] = STATE_TREE;
+ } else {
+ _forward[u] = true;
+ _pi[u] = -ART_COST;
+ _pred[u] = f;
+ _source[f] = u;
+ _target[f] = _root;
+ _cap[f] = INF;
+ _flow[f] = _supply[u];
+ _state[f] = STATE_TREE;
+ _cost[f] = ART_COST;
+ _source[e] = _root;
+ _target[e] = u;
+ _cap[e] = INF;
+ _flow[e] = 0;
+ _cost[e] = 0;
+ _state[e] = STATE_LOWER;
+ ++f;
+ }
+ }
+ _all_arc_num = f;
}
@@ -1375,18 +1438,6 @@
// Check feasibility
- if (_sum_supply < 0) {
- for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
- if (_supply[u] >= 0 && _flow[e] != 0) return INFEASIBLE;
- }
- }
- else if (_sum_supply > 0) {
- for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
- if (_supply[u] <= 0 && _flow[e] != 0) return INFEASIBLE;
- }
- }
- else {
- for (int u = 0, e = _arc_num; u != _node_num; ++u, ++e) {
- if (_flow[e] != 0) return INFEASIBLE;
- }
+ for (int e = _search_arc_num; e != _all_arc_num; ++e) {
+ if (_flow[e] != 0) return INFEASIBLE;
}
@@ -1402,4 +1453,28 @@
}
}
+
+ // Shift potentials to meet the requirements of the GEQ/LEQ type
+ // optimality conditions
+ if (_sum_supply == 0) {
+ if (_stype == GEQ) {
+ Cost max_pot = std::numeric_limits::min();
+ for (int i = 0; i != _node_num; ++i) {
+ if (_pi[i] > max_pot) max_pot = _pi[i];
+ }
+ if (max_pot > 0) {
+ for (int i = 0; i != _node_num; ++i)
+ _pi[i] -= max_pot;
+ }
+ } else {
+ Cost min_pot = std::numeric_limits::max();
+ for (int i = 0; i != _node_num; ++i) {
+ if (_pi[i] < min_pot) min_pot = _pi[i];
+ }
+ if (min_pot < 0) {
+ for (int i = 0; i != _node_num; ++i)
+ _pi[i] -= min_pot;
+ }
+ }
+ }
return OPTIMAL;