[Lemon-commits] Alpar Juttner: Merge #340
Lemon HG
hg at lemon.cs.elte.hu
Fri Feb 26 14:20:22 CET 2010
details: http://lemon.cs.elte.hu/hg/lemon/rev/2914b6f0fde0
changeset: 911:2914b6f0fde0
user: Alpar Juttner <alpar [at] cs.elte.hu>
date: Fri Feb 26 14:00:20 2010 +0100
description:
Merge #340
diffstat:
lemon/capacity_scaling.h | 15 +-
lemon/cost_scaling.h | 376 ++++++++++++++++++++++++++++++----------------
lemon/cycle_canceling.h | 9 +-
lemon/network_simplex.h | 138 ++++++++++++++--
4 files changed, 377 insertions(+), 161 deletions(-)
diffs (truncated from 915 to 300 lines):
diff --git a/lemon/capacity_scaling.h b/lemon/capacity_scaling.h
--- a/lemon/capacity_scaling.h
+++ b/lemon/capacity_scaling.h
@@ -139,9 +139,10 @@
TEMPLATE_DIGRAPH_TYPEDEFS(GR);
typedef std::vector<int> IntVector;
- typedef std::vector<char> BoolVector;
typedef std::vector<Value> ValueVector;
typedef std::vector<Cost> CostVector;
+ typedef std::vector<char> BoolVector;
+ // Note: vector<char> is used instead of vector<bool> for efficiency reasons
private:
@@ -798,15 +799,15 @@
// Initialize delta value
if (_factor > 1) {
// With scaling
- Value max_sup = 0, max_dem = 0;
- for (int i = 0; i != _node_num; ++i) {
+ Value max_sup = 0, max_dem = 0, max_cap = 0;
+ for (int i = 0; i != _root; ++i) {
Value ex = _excess[i];
if ( ex > max_sup) max_sup = ex;
if (-ex > max_dem) max_dem = -ex;
- }
- Value max_cap = 0;
- for (int j = 0; j != _res_arc_num; ++j) {
- if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
+ int last_out = _first_out[i+1] - 1;
+ for (int j = _first_out[i]; j != last_out; ++j) {
+ if (_res_cap[j] > max_cap) max_cap = _res_cap[j];
+ }
}
max_sup = std::min(std::min(max_sup, max_dem), max_cap);
for (_delta = 1; 2 * _delta <= max_sup; _delta *= 2) ;
diff --git a/lemon/cost_scaling.h b/lemon/cost_scaling.h
--- a/lemon/cost_scaling.h
+++ b/lemon/cost_scaling.h
@@ -201,10 +201,11 @@
TEMPLATE_DIGRAPH_TYPEDEFS(GR);
typedef std::vector<int> IntVector;
- typedef std::vector<char> BoolVector;
typedef std::vector<Value> ValueVector;
typedef std::vector<Cost> CostVector;
typedef std::vector<LargeCost> LargeCostVector;
+ typedef std::vector<char> BoolVector;
+ // Note: vector<char> is used instead of vector<bool> for efficiency reasons
private:
@@ -248,6 +249,7 @@
// Parameters of the problem
bool _have_lower;
Value _sum_supply;
+ int _sup_node_num;
// Data structures for storing the digraph
IntNodeMap _node_id;
@@ -276,6 +278,12 @@
LargeCost _epsilon;
int _alpha;
+ IntVector _buckets;
+ IntVector _bucket_next;
+ IntVector _bucket_prev;
+ IntVector _rank;
+ int _max_rank;
+
// Data for a StaticDigraph structure
typedef std::pair<int, int> IntPair;
StaticDigraph _sgr;
@@ -828,6 +836,11 @@
}
}
+ _sup_node_num = 0;
+ for (NodeIt n(_graph); n != INVALID; ++n) {
+ if (sup[n] > 0) ++_sup_node_num;
+ }
+
// Find a feasible flow using Circulation
Circulation<Digraph, ConstMap<Arc, Value>, ValueArcMap, ValueNodeMap>
circ(_graph, low, cap, sup);
@@ -862,7 +875,7 @@
}
for (int a = _first_out[_root]; a != _res_arc_num; ++a) {
int ra = _reverse[a];
- _res_cap[a] = 1;
+ _res_cap[a] = 0;
_res_cap[ra] = 0;
_cost[a] = 0;
_cost[ra] = 0;
@@ -876,7 +889,14 @@
void start(Method method) {
// Maximum path length for partial augment
const int MAX_PATH_LENGTH = 4;
-
+
+ // Initialize data structures for buckets
+ _max_rank = _alpha * _res_node_num;
+ _buckets.resize(_max_rank);
+ _bucket_next.resize(_res_node_num + 1);
+ _bucket_prev.resize(_res_node_num + 1);
+ _rank.resize(_res_node_num + 1);
+
// Execute the algorithm
switch (method) {
case PUSH:
@@ -915,63 +935,175 @@
}
}
}
+
+ // Initialize a cost scaling phase
+ void initPhase() {
+ // Saturate arcs not satisfying the optimality condition
+ for (int u = 0; u != _res_node_num; ++u) {
+ int last_out = _first_out[u+1];
+ LargeCost pi_u = _pi[u];
+ for (int a = _first_out[u]; a != last_out; ++a) {
+ int v = _target[a];
+ if (_res_cap[a] > 0 && _cost[a] + pi_u - _pi[v] < 0) {
+ Value delta = _res_cap[a];
+ _excess[u] -= delta;
+ _excess[v] += delta;
+ _res_cap[a] = 0;
+ _res_cap[_reverse[a]] += delta;
+ }
+ }
+ }
+
+ // Find active nodes (i.e. nodes with positive excess)
+ for (int u = 0; u != _res_node_num; ++u) {
+ if (_excess[u] > 0) _active_nodes.push_back(u);
+ }
+
+ // Initialize the next arcs
+ for (int u = 0; u != _res_node_num; ++u) {
+ _next_out[u] = _first_out[u];
+ }
+ }
+
+ // Early termination heuristic
+ bool earlyTermination() {
+ const double EARLY_TERM_FACTOR = 3.0;
+
+ // Build a static residual graph
+ _arc_vec.clear();
+ _cost_vec.clear();
+ for (int j = 0; j != _res_arc_num; ++j) {
+ if (_res_cap[j] > 0) {
+ _arc_vec.push_back(IntPair(_source[j], _target[j]));
+ _cost_vec.push_back(_cost[j] + 1);
+ }
+ }
+ _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
+
+ // Run Bellman-Ford algorithm to check if the current flow is optimal
+ BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
+ bf.init(0);
+ bool done = false;
+ int K = int(EARLY_TERM_FACTOR * std::sqrt(double(_res_node_num)));
+ for (int i = 0; i < K && !done; ++i) {
+ done = bf.processNextWeakRound();
+ }
+ return done;
+ }
+
+ // Global potential update heuristic
+ void globalUpdate() {
+ int bucket_end = _root + 1;
+
+ // Initialize buckets
+ for (int r = 0; r != _max_rank; ++r) {
+ _buckets[r] = bucket_end;
+ }
+ Value total_excess = 0;
+ for (int i = 0; i != _res_node_num; ++i) {
+ if (_excess[i] < 0) {
+ _rank[i] = 0;
+ _bucket_next[i] = _buckets[0];
+ _bucket_prev[_buckets[0]] = i;
+ _buckets[0] = i;
+ } else {
+ total_excess += _excess[i];
+ _rank[i] = _max_rank;
+ }
+ }
+ if (total_excess == 0) return;
+
+ // Search the buckets
+ int r = 0;
+ for ( ; r != _max_rank; ++r) {
+ while (_buckets[r] != bucket_end) {
+ // Remove the first node from the current bucket
+ int u = _buckets[r];
+ _buckets[r] = _bucket_next[u];
+
+ // Search the incomming arcs of u
+ LargeCost pi_u = _pi[u];
+ int last_out = _first_out[u+1];
+ for (int a = _first_out[u]; a != last_out; ++a) {
+ int ra = _reverse[a];
+ if (_res_cap[ra] > 0) {
+ int v = _source[ra];
+ int old_rank_v = _rank[v];
+ if (r < old_rank_v) {
+ // Compute the new rank of v
+ LargeCost nrc = (_cost[ra] + _pi[v] - pi_u) / _epsilon;
+ int new_rank_v = old_rank_v;
+ if (nrc < LargeCost(_max_rank))
+ new_rank_v = r + 1 + int(nrc);
+
+ // Change the rank of v
+ if (new_rank_v < old_rank_v) {
+ _rank[v] = new_rank_v;
+ _next_out[v] = _first_out[v];
+
+ // Remove v from its old bucket
+ if (old_rank_v < _max_rank) {
+ if (_buckets[old_rank_v] == v) {
+ _buckets[old_rank_v] = _bucket_next[v];
+ } else {
+ _bucket_next[_bucket_prev[v]] = _bucket_next[v];
+ _bucket_prev[_bucket_next[v]] = _bucket_prev[v];
+ }
+ }
+
+ // Insert v to its new bucket
+ _bucket_next[v] = _buckets[new_rank_v];
+ _bucket_prev[_buckets[new_rank_v]] = v;
+ _buckets[new_rank_v] = v;
+ }
+ }
+ }
+ }
+
+ // Finish search if there are no more active nodes
+ if (_excess[u] > 0) {
+ total_excess -= _excess[u];
+ if (total_excess <= 0) break;
+ }
+ }
+ if (total_excess <= 0) break;
+ }
+
+ // Relabel nodes
+ for (int u = 0; u != _res_node_num; ++u) {
+ int k = std::min(_rank[u], r);
+ if (k > 0) {
+ _pi[u] -= _epsilon * k;
+ _next_out[u] = _first_out[u];
+ }
+ }
+ }
/// Execute the algorithm performing augment and relabel operations
void startAugment(int max_length = std::numeric_limits<int>::max()) {
// Paramters for heuristics
- const int BF_HEURISTIC_EPSILON_BOUND = 1000;
- const int BF_HEURISTIC_BOUND_FACTOR = 3;
+ const int EARLY_TERM_EPSILON_LIMIT = 1000;
+ const double GLOBAL_UPDATE_FACTOR = 3.0;
+ const int global_update_freq = int(GLOBAL_UPDATE_FACTOR *
+ (_res_node_num + _sup_node_num * _sup_node_num));
+ int next_update_limit = global_update_freq;
+
+ int relabel_cnt = 0;
+
// Perform cost scaling phases
- IntVector pred_arc(_res_node_num);
- std::vector<int> path_nodes;
+ std::vector<int> path;
for ( ; _epsilon >= 1; _epsilon = _epsilon < _alpha && _epsilon > 1 ?
1 : _epsilon / _alpha )
{
- // "Early Termination" heuristic: use Bellman-Ford algorithm
- // to check if the current flow is optimal
- if (_epsilon <= BF_HEURISTIC_EPSILON_BOUND) {
- _arc_vec.clear();
- _cost_vec.clear();
- for (int j = 0; j != _res_arc_num; ++j) {
- if (_res_cap[j] > 0) {
- _arc_vec.push_back(IntPair(_source[j], _target[j]));
- _cost_vec.push_back(_cost[j] + 1);
- }
- }
- _sgr.build(_res_node_num, _arc_vec.begin(), _arc_vec.end());
-
- BellmanFord<StaticDigraph, LargeCostArcMap> bf(_sgr, _cost_map);
- bf.init(0);
- bool done = false;
- int K = int(BF_HEURISTIC_BOUND_FACTOR * sqrt(_res_node_num));
- for (int i = 0; i < K && !done; ++i)
- done = bf.processNextWeakRound();
- if (done) break;
More information about the Lemon-commits
mailing list