From 0d6ffe6b31faa004f752f47070c414fa03b34f9e Mon Sep 17 00:00:00 2001 From: Anh-Dung Phan Date: Mon, 11 Nov 2013 08:51:52 +0100 Subject: [PATCH] Implement three pivot rules --- src/smt/network_flow.h | 195 ++++++++++++++++++++++++++++++++++-- src/smt/network_flow_def.h | 45 ++++----- src/smt/spanning_tree_def.h | 34 ++++--- 3 files changed, 226 insertions(+), 48 deletions(-) diff --git a/src/smt/network_flow.h b/src/smt/network_flow.h index 44861575a..5390f53ea 100644 --- a/src/smt/network_flow.h +++ b/src/smt/network_flow.h @@ -23,8 +23,6 @@ Notes: It remains unclear how to convert DL assignment to a basic feasible solution of Network Simplex. A naive approach is to run an algorithm on max flow in order to get a spanning tree. - - The network_simplex class hasn't had multiple pivoting strategies yet. --*/ #ifndef _NETWORK_FLOW_H_ @@ -36,20 +34,207 @@ Notes: namespace smt { + enum pivot_rule { + // First eligible edge pivot rule + // Edges are traversed in a wraparound fashion + FIRST_ELIGIBLE, + // Best eligible edge pivot rule + // The best edge is selected in every iteration + BEST_ELIGIBLE, + // Candidate list pivot rule + // Major iterations: candidate list is built from eligible edges (in a wraparound way) + // Minor iterations: the best edge is selected from the list + CANDIDATE_LIST + }; + // Solve minimum cost flow problem using Network Simplex algorithm template class network_flow : private Ext { + private: enum edge_state { LOWER = 1, BASIS = 0, UPPER = -1 }; + typedef dl_var node; typedef dl_edge edge; typedef dl_graph graph; typedef typename Ext::numeral numeral; typedef typename Ext::fin_numeral fin_numeral; + class pivot_rule_impl { + protected: + graph & m_graph; + svector & m_states; + vector & m_potentials; + edge_id & m_enter_id; + + public: + pivot_rule_impl() {} + pivot_rule_impl(graph & g, vector & potentials, + svector & states, edge_id & enter_id) + : m_graph(g), + m_potentials(potentials), + m_states(states), + m_enter_id(enter_id) { + } + bool choose_entering_edge() {return false;}; + }; + + class first_eligible_pivot : pivot_rule_impl { + private: + edge_id m_next_edge; + + public: + first_eligible_pivot(graph & g, vector & potentials, + svector & states, edge_id & enter_id) : + pivot_rule_impl(g, potentials, states, enter_id), + m_next_edge(0) { + } + + bool choose_entering_edge() { + TRACE("network_flow", tout << "choose_entering_edge...\n";); + unsigned num_edges = m_graph.get_num_edges(); + for (unsigned i = m_next_edge; i < m_next_edge + num_edges; ++i) { + edge_id id = (i >= num_edges) ? (i - num_edges) : i; + node src = m_graph.get_source(id); + node tgt = m_graph.get_target(id); + if (m_states[id] != BASIS) { + numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id); + if (cost.is_pos()) { + m_enter_id = id; + TRACE("network_flow", { + tout << "Found entering edge " << id << " between node "; + tout << src << " and node " << tgt << " with reduced cost = " << cost << "...\n"; + }); + return true; + } + } + } + TRACE("network_flow", tout << "Found no entering edge...\n";); + return false; + }; + }; + + class best_eligible_pivot : pivot_rule_impl { + public: + best_eligible_pivot(graph & g, vector & potentials, + svector & states, edge_id & enter_id) : + pivot_rule_impl(g, potentials, states, enter_id) { + } + + bool choose_entering_edge() { + TRACE("network_flow", tout << "choose_entering_edge...\n";); + unsigned num_edges = m_graph.get_num_edges(); + numeral max = numeral::zero(); + for (unsigned i = 0; i < num_edges; ++i) { + node src = m_graph.get_source(i); + node tgt = m_graph.get_target(i); + if (m_states[i] != BASIS) { + numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(i); + if (cost > max) { + max = cost; + m_enter_id = i; + } + } + } + if (max.is_pos()) { + TRACE("network_flow", { + tout << "Found entering edge " << m_enter_id << " between node "; + tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id); + tout << " with reduced cost = " << max << "...\n"; + }); + return true; + } + TRACE("network_flow", tout << "Found no entering edge...\n";); + return false; + }; + }; + + class candidate_list_pivot : pivot_rule_impl { + private: + edge_id m_next_edge; + svector m_candidates; + unsigned num_candidates; + static const unsigned NUM_CANDIDATES = 10; + + public: + candidate_list_pivot(graph & g, vector & potentials, + svector & states, edge_id & enter_id) : + pivot_rule_impl(g, potentials, states, enter_id), + m_next_edge(0), + num_candidates(NUM_CANDIDATES), + m_candidates(num_candidates) { + } + + bool choose_entering_edge() { + if (m_candidates.empty()) { + // Build the candidate list + unsigned num_edges = m_graph.get_num_edges(); + numeral max = numeral::zero(); + unsigned count = 0; + for (unsigned i = m_next_edge; i < m_next_edge + num_edges; ++i) { + edge_id id = (i >= num_edges) ? i - num_edges : i; + node src = m_graph.get_source(id); + node tgt = m_graph.get_target(id); + if (m_states[id] != BASIS) { + numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id); + if (cost.is_pos()) { + m_candidates[count++] = id; + if (cost > max) { + max = cost; + m_enter_id = id; + } + } + if (count >= num_candidates) break; + } + } + m_next_edge = m_enter_id; + if (max.is_pos()) { + TRACE("network_flow", { + tout << "Found entering edge " << m_enter_id << " between node "; + tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id); + tout << " with reduced cost = " << max << "...\n"; + }); + return true; + } + TRACE("network_flow", tout << "Found no entering edge...\n";); + return false; + } + else { + numeral max = numeral::zero(); + unsigned last = m_candidates.size(); + for (unsigned i = 0; i < last; ++i) { + edge_id id = m_candidates[i]; + node src = m_graph.get_source(id); + node tgt = m_graph.get_target(id); + if (m_states[id] != BASIS) { + numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id); + if (cost > max) { + max = cost; + m_enter_id = id; + } + // Remove stale candidates + if (!cost.is_pos()) { + m_candidates[i] = m_candidates[--last]; + } + } + } + if (max.is_pos()) { + TRACE("network_flow", { + tout << "Found entering edge " << m_enter_id << " between node "; + tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id); + tout << " with reduced cost = " << max << "...\n"; + }); + return true; + } + TRACE("network_flow", tout << "Found no entering edge...\n";); + return false; + } + }; + }; + graph m_graph; thread_spanning_tree m_tree; @@ -76,9 +261,7 @@ namespace smt { void update_flows(); - // If all reduced costs are non-negative, return false since the current spanning tree is optimal - // Otherwise return true and update m_entering_edge - bool choose_entering_edge(); + bool choose_entering_edge(pivot_rule pr); // Send as much flow as possible around the cycle, the first basic edge with flow 0 will leave // Return false if the problem is unbounded @@ -99,7 +282,7 @@ namespace smt { // Minimize cost flows // Return true if found an optimal solution, and return false if unbounded - bool min_cost(); + bool min_cost(pivot_rule pr = FIRST_ELIGIBLE); // Compute the optimal solution numeral get_optimal_solution(vector & result, bool is_dual); diff --git a/src/smt/network_flow_def.h b/src/smt/network_flow_def.h index d31289dfc..515d49933 100644 --- a/src/smt/network_flow_def.h +++ b/src/smt/network_flow_def.h @@ -129,30 +129,6 @@ namespace smt { TRACE("network_flow", tout << pp_vector("Flows", m_flows, true);); } - template - bool network_flow::choose_entering_edge() { - TRACE("network_flow", tout << "choose_entering_edge...\n";); - unsigned num_edges = m_graph.get_num_edges(); - for (unsigned i = 0; i < num_edges; ++i) { - node src = m_graph.get_source(i); - node tgt = m_graph.get_target(i); - if (m_states[i] != BASIS) { - numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(i); - // TODO: add multiple pivoting strategies - if (cost.is_pos()) { - m_enter_id = i; - TRACE("network_flow", { - tout << "Found entering edge " << i << " between node "; - tout << src << " and node " << tgt << " with reduced cost = " << cost << "...\n"; - }); - return true; - } - } - } - TRACE("network_flow", tout << "Found no entering edge...\n";); - return false; - } - template bool network_flow::choose_leaving_edge() { TRACE("network_flow", tout << "choose_leaving_edge...\n";); @@ -191,12 +167,29 @@ namespace smt { m_tree.update(m_enter_id, m_leave_id); } + // FIXME: should declare pivot as a pivot_rule_impl and refactor + template + bool network_flow::choose_entering_edge(pivot_rule pr) { + if (pr == FIRST_ELIGIBLE) { + first_eligible_pivot pivot(m_graph, m_potentials, m_states, m_enter_id); + return pivot.choose_entering_edge(); + } + else if (pr == BEST_ELIGIBLE) { + best_eligible_pivot pivot(m_graph, m_potentials, m_states, m_enter_id); + return pivot.choose_entering_edge(); + } + else { + candidate_list_pivot pivot(m_graph, m_potentials, m_states, m_enter_id); + return pivot.choose_entering_edge(); + } + } + // Minimize cost flows // Return true if found an optimal solution, and return false if unbounded template - bool network_flow::min_cost() { + bool network_flow::min_cost(pivot_rule pr) { initialize(); - while (choose_entering_edge()) { + while (choose_entering_edge(pr)) { bool bounded = choose_leaving_edge(); if (!bounded) return false; update_flows(); diff --git a/src/smt/spanning_tree_def.h b/src/smt/spanning_tree_def.h index 2be593842..be668d91d 100644 --- a/src/smt/spanning_tree_def.h +++ b/src/smt/spanning_tree_def.h @@ -162,7 +162,7 @@ namespace smt { tout << u << ", " << v << ") leaves\n"; }); - node old_pred = m_pred[q]; + node old_pred = m_pred[q]; // Update stem nodes from q to v if (q != v) { for (node n = q; n != u; ) { @@ -175,18 +175,23 @@ namespace smt { old_pred = next_old_pred; } } - else { - node x = get_final(p); - node y = m_thread[x]; - node z = get_final(q); - node t = m_thread[get_final(v)]; - node r = find_rev_thread(v); - m_thread[z] = y; - m_thread[x] = q; - m_thread[r] = t; - } - m_pred[q] = p; + // Old threads: alpha -> q -*-> f(q) -> beta | p -*-> f(p) -> gamma + // New threads: alpha -> beta | p -*-> f(p) -> q -*-> f(q) -> gamma + + node f_p = get_final(p); + node f_q = get_final(q); + node alpha = find_rev_thread(q); + node beta = m_thread[f_q]; + node gamma = m_thread[f_p]; + + if (q != gamma) { + m_thread[alpha] = beta; + m_thread[f_p] = q; + m_thread[f_q] = gamma; + } + + m_pred[q] = p; m_tree[q] = enter_id; m_root_t2 = q; @@ -211,7 +216,6 @@ namespace smt { Spanning tree of m_graph + root is represented using: svector m_states; edge_id |-> edge_state - svector m_upwards; node |-> bool svector m_pred; node |-> node svector m_depth; node |-> int svector m_thread; node |-> node @@ -224,9 +228,7 @@ namespace smt { m_thread is a linked list traversing all nodes. Furthermore, the nodes linked in m_thread follows a depth-first traversal order. - - m_upwards direction of edge from i to m_pred[i] m_graph - + */ template bool thread_spanning_tree::check_well_formed() {