mirrors/z3
3
0
Fork 0
mirror of https://github.com/Z3Prover/z3 synced 2025-06-18 11:58:31 +00:00
Code Activity

Create a separate class for spanning tree

Browse source 
Remarks:

1. Templates should be in header files only

2. Should pass in svector<_> instead of returning a local one
This commit is contained in:
Anh-Dung Phan 2013-11-06 17:42:09 -08:00
parent 034b33b6da
commit f7fdf134fd
5 changed files with 586 additions and 510 deletions
Download patch file Download diff file Expand all files Collapse all files

36
src/smt/network_flow.h
View file

@ -32,12 +32,10 @@ Notes:
#include"inf_rational.h" #include"inf_rational.h"
#include"diff_logic.h" #include"diff_logic.h"
#include"spanning_tree.h"
namespace smt { namespace smt {
template<typename TV>
std::string pp_vector(std::string const & label, TV v, bool has_header = false);
// Solve minimum cost flow problem using Network Simplex algorithm // Solve minimum cost flow problem using Network Simplex algorithm
template<typename Ext> template<typename Ext>
class network_flow : private Ext { class network_flow : private Ext {
@ -51,7 +49,9 @@ namespace smt {
typedef dl_graph<Ext> graph; typedef dl_graph<Ext> graph;
typedef typename Ext::numeral numeral; typedef typename Ext::numeral numeral;
typedef typename Ext::fin_numeral fin_numeral; typedef typename Ext::fin_numeral fin_numeral;
graph m_graph; graph m_graph;
thread_spanning_tree tree;
// Denote supply/demand b_i on node i // Denote supply/demand b_i on node i
vector<fin_numeral> m_balances; vector<fin_numeral> m_balances;
@ -67,18 +67,7 @@ namespace smt {
svector<edge_state> m_states; svector<edge_state> m_states;
// m_upwards[i] is true if the corresponding edge unsigned m_step;
// (i, m_pred[i]) points upwards (pointing toward the root node)
svector<bool> m_upwards;
// Store the parent of a node i in the spanning tree
svector<node> m_pred;
// Store the number of edge on the path from node i to the root
svector<int> m_depth;
// Store the pointer from node i to the next node in depth-first search order
svector<node> m_thread;
// Store a final node of the sub tree rooted at node i
svector<node> m_final;
edge_id m_entering_edge; edge_id m_entering_edge;
edge_id m_leaving_edge; edge_id m_leaving_edge;
@ -86,14 +75,11 @@ namespace smt {
optional<numeral> m_delta; optional<numeral> m_delta;
bool m_in_edge_dir; bool m_in_edge_dir;
unsigned m_step;
// Initialize the network with a feasible spanning tree // Initialize the network with a feasible spanning tree
void initialize(); void initialize();
edge_id get_edge_id(dl_var source, dl_var target) const; edge_id get_edge_id(dl_var source, dl_var target) const;
void update_potentials(); void update_potentials();
void update_flows(); void update_flows();
@ -113,22 +99,8 @@ namespace smt {
bool edge_in_tree(edge_id id) const; bool edge_in_tree(edge_id id) const;
bool edge_in_tree(node src, node dst) const; bool edge_in_tree(node src, node dst) const;
bool is_ancestor_of(node ancestor, node child) const;
/**
\brief find node that points to 'n' in m_thread
*/
node find_rev_thread(node n) const;
void fix_depth(node start, node end);
void swap_order(node q, node v);
bool check_well_formed(); bool check_well_formed();
bool is_preorder_traversal(node start, node end);
node get_final(int start);
public: public:
network_flow(graph & g, vector<fin_numeral> const & balances); network_flow(graph & g, vector<fin_numeral> const & balances);

560
src/smt/network_flow_def.h
View file

@ -25,24 +25,6 @@ Notes:
namespace smt { namespace smt {
template<typename TV>
std::string pp_vector(std::string const & label, TV v, bool has_header) {
std::ostringstream oss;
if (has_header) {
oss << "Index ";
for (unsigned i = 0; i < v.size(); ++i) {
oss << i << " ";
}
oss << std::endl;
}
oss << label << " ";
for (unsigned i = 0; i < v.size(); ++i) {
oss << v[i] << " ";
}
oss << std::endl;
return oss.str();
}
template<typename Ext> template<typename Ext>
network_flow<Ext>::network_flow(graph & g, vector<fin_numeral> const & balances) : network_flow<Ext>::network_flow(graph & g, vector<fin_numeral> const & balances) :
m_balances(balances) { m_balances(balances) {
@ -65,12 +47,7 @@ namespace smt {
m_balances.resize(num_nodes); m_balances.resize(num_nodes);
m_potentials.resize(num_nodes); m_potentials.resize(num_nodes);
m_upwards.resize(num_nodes); tree = thread_spanning_tree();
m_pred.resize(num_nodes);
m_depth.resize(num_nodes);
m_thread.resize(num_nodes);
m_step = 0;
} }
template<typename Ext> template<typename Ext>
@ -79,13 +56,10 @@ namespace smt {
// Create an artificial root node to construct initial spanning tree // Create an artificial root node to construct initial spanning tree
unsigned num_nodes = m_graph.get_num_nodes(); unsigned num_nodes = m_graph.get_num_nodes();
unsigned num_edges = m_graph.get_num_edges(); unsigned num_edges = m_graph.get_num_edges();
node root = num_nodes;
m_pred[root] = -1;
m_depth[root] = 0;
m_thread[root] = 0;
m_potentials[root] = numeral::zero();
node root = num_nodes;
m_graph.init_var(root); m_graph.init_var(root);
m_potentials[root] = numeral::zero();
fin_numeral sum_supply = fin_numeral::zero(); fin_numeral sum_supply = fin_numeral::zero();
for (unsigned i = 0; i < m_balances.size(); ++i) { for (unsigned i = 0; i < m_balances.size(); ++i) {
@ -99,31 +73,31 @@ namespace smt {
m_states.fill(LOWER); m_states.fill(LOWER);
// Create artificial edges from/to root node to/from other nodes and initialize the spanning tree // Create artificial edges from/to root node to/from other nodes and initialize the spanning tree
svector<bool> upwards(num_nodes, false);
for (unsigned i = 0; i < num_nodes; ++i) { for (unsigned i = 0; i < num_nodes; ++i) {
m_upwards[i] = !m_balances[i].is_neg(); upwards[i] = !m_balances[i].is_neg();
m_pred[i] = root;
m_depth[i] = 1;
m_thread[i] = i + 1;
m_states[num_edges + i] = BASIS; m_states[num_edges + i] = BASIS;
node src = m_upwards[i] ? i : root; node src = upwards[i] ? i : root;
node tgt = m_upwards[i] ? root : i; node tgt = upwards[i] ? root : i;
m_flows[num_edges + i] = m_upwards[i] ? m_balances[i] : -m_balances[i]; m_flows[num_edges + i] = upwards[i] ? m_balances[i] : -m_balances[i];
m_graph.add_edge(src, tgt, numeral::one(), explanation()); m_graph.add_edge(src, tgt, numeral::one(), explanation());
} }
tree.initialize(upwards, num_nodes);
// Compute initial potentials // Compute initial potentials
node u = m_thread[root]; svector<node> descendants;
while (u != root) { tree.get_descendants(root, descendants);
node v = m_pred[u]; // Skip root node
for (unsigned i = 1; i < descendants.size(); ++i) {
node u = descendants[i];
node v = tree.get_parent(u);
edge_id e_id = get_edge_id(u, v); edge_id e_id = get_edge_id(u, v);
m_potentials[u] = m_potentials[v] + (m_upwards[u] ? - m_graph.get_weight(e_id) : m_graph.get_weight(e_id)); m_potentials[u] = m_potentials[v] + (tree.get_arc_direction(u) ? - m_graph.get_weight(e_id) : m_graph.get_weight(e_id));
u = m_thread[u];
} }
TRACE("network_flow", { TRACE("network_flow", {
tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread); tout << pp_vector("Potentials", m_potentials, true) << pp_vector("Flows", m_flows);
tout << pp_vector("Depths", m_depth) << pp_vector("Upwards", m_upwards);
tout << pp_vector("Potentials", m_potentials) << pp_vector("Flows", m_flows);
}); });
TRACE("network_flow", tout << "Spanning tree:\n" << display_spanning_tree();); TRACE("network_flow", tout << "Spanning tree:\n" << display_spanning_tree(););
SASSERT(check_well_formed()); SASSERT(check_well_formed());
@ -131,9 +105,9 @@ namespace smt {
template<typename Ext> template<typename Ext>
edge_id network_flow<Ext>::get_edge_id(dl_var source, dl_var target) const { edge_id network_flow<Ext>::get_edge_id(dl_var source, dl_var target) const {
// m_upwards[source] decides which node is the real source // tree.get_arc_direction(source) decides which node is the real source
edge_id id; edge_id id;
VERIFY(m_upwards[source] ? m_graph.get_edge_id(source, target, id) : m_graph.get_edge_id(target, source, id)); VERIFY(tree.get_arc_direction(source) ? m_graph.get_edge_id(source, target, id) : m_graph.get_edge_id(target, source, id));
return id; return id;
} }
@ -143,9 +117,11 @@ namespace smt {
node src = m_graph.get_source(m_entering_edge); node src = m_graph.get_source(m_entering_edge);
node tgt = m_graph.get_target(m_entering_edge); node tgt = m_graph.get_target(m_entering_edge);
numeral cost = m_graph.get_weight(m_entering_edge); numeral cost = m_graph.get_weight(m_entering_edge);
numeral change = m_upwards[src] ? (-cost - m_potentials[src] + m_potentials[tgt]) : (cost + m_potentials[src] - m_potentials[tgt]); numeral change = tree.get_arc_direction(src) ? (-cost - m_potentials[src] + m_potentials[tgt]) : (cost + m_potentials[src] - m_potentials[tgt]);
node last = m_thread[get_final(src)]; svector<node> descendants;
for (node u = src; u != last; u = m_thread[u]) { tree.get_descendants(src, descendants);
for (unsigned i = 0; i < descendants.size(); ++i) {
node u = descendants[i];
m_potentials[u] += change; m_potentials[u] += change;
} }
TRACE("network_flow", tout << pp_vector("Potentials", m_potentials, true);); TRACE("network_flow", tout << pp_vector("Potentials", m_potentials, true););
@ -157,14 +133,20 @@ namespace smt {
numeral val = fin_numeral(m_states[m_entering_edge]) * (*m_delta); numeral val = fin_numeral(m_states[m_entering_edge]) * (*m_delta);
m_flows[m_entering_edge] += val; m_flows[m_entering_edge] += val;
node source = m_graph.get_source(m_entering_edge); node source = m_graph.get_source(m_entering_edge);
for (unsigned u = source; u != m_join_node; u = m_pred[u]) { svector<node> ancestors;
edge_id e_id = get_edge_id(u, m_pred[u]); tree.get_ancestors(source, ancestors);
m_flows[e_id] += m_upwards[u] ? -val : val; for (unsigned i = 0; i < ancestors.size() && ancestors[i] != m_join_node; ++i) {
node u = ancestors[i];
edge_id e_id = get_edge_id(u, tree.get_parent(u));
m_flows[e_id] += tree.get_arc_direction(u) ? -val : val;
} }
node target = m_graph.get_target(m_entering_edge); node target = m_graph.get_target(m_entering_edge);
for (unsigned u = target; u != m_join_node; u = m_pred[u]) { tree.get_ancestors(target,ancestors);
edge_id e_id = get_edge_id(u, m_pred[u]); for (unsigned i = 0; i < ancestors.size() && ancestors[i] != m_join_node; ++i) {
m_flows[e_id] += m_upwards[u] ? val : -val; node u = ancestors[i];
edge_id e_id = get_edge_id(u, tree.get_parent(u));
m_flows[e_id] += tree.get_arc_direction(u) ? val : -val;
} }
TRACE("network_flow", tout << pp_vector("Flows", m_flows, true);); TRACE("network_flow", tout << pp_vector("Flows", m_flows, true););
} }
@ -202,41 +184,37 @@ namespace smt {
if (m_states[m_entering_edge] == UPPER) { if (m_states[m_entering_edge] == UPPER) {
std::swap(source, target); std::swap(source, target);
} }
node u = source, v = target;
while (u != v) { m_join_node = tree.get_common_ancestor(source, target);
if (m_depth[u] > m_depth[v])
u = m_pred[u];
else if (m_depth[u] < m_depth[v])
v = m_pred[v];
else {
u = m_pred[u];
v = m_pred[v];
}
}
// Found first common ancestor of source and target
m_join_node = u;
TRACE("network_flow", tout << "Found join node " << m_join_node << std::endl;); TRACE("network_flow", tout << "Found join node " << m_join_node << std::endl;);
m_delta.set_invalid(); m_delta.set_invalid();
node src, tgt; node src, tgt;
// Send flows along the path from source to the ancestor // Send flows along the path from source to the ancestor
for (unsigned u = source; u != m_join_node; u = m_pred[u]) { svector<node> ancestors;
edge_id e_id = get_edge_id(u, m_pred[u]); tree.get_ancestors(source, ancestors);
if (m_upwards[u] && (!m_delta || m_flows[e_id] < *m_delta)) { for (unsigned i = 0; i < ancestors.size() && ancestors[i] != m_join_node; ++i) {
node u = ancestors[i];
edge_id e_id = get_edge_id(u, tree.get_parent(u));
if (tree.get_arc_direction(u) && (!m_delta || m_flows[e_id] < *m_delta)) {
m_delta = m_flows[e_id]; m_delta = m_flows[e_id];
src = u; src = u;
tgt = m_pred[u]; tgt = tree.get_parent(u);
SASSERT(edge_in_tree(src,tgt)); SASSERT(edge_in_tree(src,tgt));
m_in_edge_dir = true; m_in_edge_dir = true;
} }
} }
// Send flows along the path from target to the ancestor // Send flows along the path from target to the ancestor
for (unsigned u = target; u != m_join_node; u = m_pred[u]) { tree.get_ancestors(target, ancestors);
edge_id e_id = get_edge_id(u, m_pred[u]); for (unsigned i = 0; i < ancestors.size() && ancestors[i] != m_join_node; ++i) {
if (!m_upwards[u] && (!m_delta || m_flows[e_id] <= *m_delta)) { node u = ancestors[i];
edge_id e_id = get_edge_id(u, tree.get_parent(u));
if (!tree.get_arc_direction(u) && (!m_delta || m_flows[e_id] <= *m_delta)) {
m_delta = m_flows[e_id]; m_delta = m_flows[e_id];
src = u; src = u;
tgt = m_pred[u]; tgt = tree.get_parent(u);
SASSERT(edge_in_tree(src,tgt)); SASSERT(edge_in_tree(src,tgt));
m_in_edge_dir = false; m_in_edge_dir = false;
} }
@ -254,196 +232,14 @@ namespace smt {
return false; return false;
} }
template<typename Ext>
bool network_flow<Ext>::is_ancestor_of(node ancestor, node child) const {
for (node n = child; n != -1; n = m_pred[n]) {
if (n == ancestor) {
return true;
}
}
return false;
}
template<typename Ext>
dl_var network_flow<Ext>::find_rev_thread(node n) const {
node ancestor = m_pred[n];
SASSERT(ancestor != -1);
while (m_thread[ancestor] != n) {
ancestor = m_thread[ancestor];
}
return ancestor;
}
template<typename Ext>
void network_flow<Ext>::fix_depth(node start, node end) {
SASSERT(m_pred[start] != -1);
m_depth[start] = m_depth[m_pred[start]]+1;
while (start != end) {
start = m_thread[start];
m_depth[start] = m_depth[m_pred[start]]+1;
}
}
/**
\brief add entering_edge, remove leaving_edge from spanning tree.
Old tree: New tree:
root root
/ \ / \
x y x y
/ \ / \ / \ / \
u s u s
| / /
v w v w
/ \ \ / \ \
z p z p
\ \ /
q q
*/
template<typename Ext> template<typename Ext>
void network_flow<Ext>::update_spanning_tree() { void network_flow<Ext>::update_spanning_tree() {
node p = m_graph.get_source(m_entering_edge); node p = m_graph.get_source(m_entering_edge);
node q = m_graph.get_target(m_entering_edge); node q = m_graph.get_target(m_entering_edge);
node u = m_graph.get_source(m_leaving_edge); node u = m_graph.get_source(m_leaving_edge);
node v = m_graph.get_target(m_leaving_edge); node v = m_graph.get_target(m_leaving_edge);
bool q_upwards = false;
// v is parent of u so T_u does not contain root node tree.update(p, q, u, v);
if (m_pred[u] == v) {
std::swap(u, v);
}
SASSERT(m_pred[v] == u);
if (is_ancestor_of(v, p)) {
std::swap(p, q);
q_upwards = true;
}
SASSERT(is_ancestor_of(v, q));
TRACE("network_flow", {
tout << "update_spanning_tree: (" << p << ", " << q << ") enters, (";
tout << u << ", " << v << ") leaves\n";
});
// Update m_pred (for nodes in the stem from q to v)
// Note: m_pred[v] == u
// Initialize m_upwards[q] = q_upwards
bool prev_upwards = q_upwards;
#if 0
for (node n = q, prev = p; n != u; ) {
SASSERT(m_pred[n] != u || n == v); // the last processed node is v
node next = m_pred[n];
m_pred[n] = prev;
std::swap(m_upwards[n], prev_upwards);
prev_upwards = !prev_upwards; // flip previous version of upwards.
prev = n;
n = next;
SASSERT(n != -1);
}
#else
node old_pred = m_pred[q];
if (q != v) {
for (node n = q; n != u; ) {
SASSERT(old_pred != u || n == v); // the last processed node is v
TRACE("network_flow", {
tout << pp_vector("Predecessors", m_pred, true);
});
SASSERT(-1 != m_pred[old_pred]);
node next_old_pred = m_pred[old_pred];
swap_order(n, old_pred);
std::swap(m_upwards[n], prev_upwards);
prev_upwards = !prev_upwards; // flip previous version of upwards.
n = old_pred;
old_pred = next_old_pred;
}
}
m_pred[q] = p;
#endif
// m_thread were updated.
// update the depth.
fix_depth(q, get_final(q));
TRACE("network_flow", {
tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread);
tout << pp_vector("Depths", m_depth) << pp_vector("Upwards", m_upwards);
});
}
/**
swap v and q in tree.
- fixup m_thread
- fixup m_pred
Case 1: final(q) == final(v)
-------
Old thread: prev -> v -*-> alpha -> q -*-> final(q) -> next
New thread: prev -> q -*-> final(q) -> v -*-> alpha -> next
Case 2: final(q) != final(v)
-------
Old thread: prev -> v -*-> alpha -> q -*-> final(q) -> beta -*-> final(v) -> next
New thread: prev -> q -*-> final(q) -> v -*-> alpha -> beta -*-> final(v) -> next
*/
template<typename Ext>
void network_flow<Ext>::swap_order(node q, node v) {
SASSERT(q != v);
SASSERT(m_pred[q] == v);
SASSERT(is_preorder_traversal(v, get_final(v)));
node prev = find_rev_thread(v);
node final_q = get_final(q);
node final_v = get_final(v);
node next = m_thread[final_v];
node alpha = find_rev_thread(q);
if (final_q == final_v) {
m_thread[final_q] = v;
m_thread[alpha] = next;
}
else {
node beta = m_thread[final_q];
m_thread[final_q] = v;
m_thread[alpha] = beta;
}
m_thread[prev] = q;
m_pred[v] = q;
SASSERT(is_preorder_traversal(q, get_final(q)));
}
template<typename Ext>
std::string network_flow<Ext>::display_spanning_tree() {
++m_step;;
std::ostringstream oss;
std::string prefix = "T";
prefix.append(std::to_string(m_step));
prefix.append("_");
unsigned root = m_graph.get_num_nodes() - 1;
for (unsigned i = 0; i < root; ++i) {
oss << prefix << i << "[shape=circle,label=\"" << prefix << i << " [";
oss << m_potentials[i] << "/" << m_balances[i] << "]\"];\n";
}
oss << prefix << root << "[shape=doublecircle,label=\"" << prefix << root << " [";
oss << m_potentials[root] << "/" << m_balances[root] << "]\"];\n";
vector<edge> const & es = m_graph.get_all_edges();
for (unsigned i = 0; i < es.size(); ++i) {
edge const & e = es[i];
oss << prefix << e.get_source() << " -> " << prefix << e.get_target();
if (m_states[i] == BASIS) {
oss << "[color=red,penwidth=3.0,label=\"" << m_flows[i] << "/" << e.get_weight() << "\"];\n";
}
else {
oss << "[label=\"" << m_flows[i] << "/" << e.get_weight() << "\"];\n";
}
}
oss << std::endl;
return oss.str();
} }
// Minimize cost flows // Minimize cost flows
@ -494,42 +290,6 @@ namespace smt {
return m_objective_value; return m_objective_value;
} }
static unsigned find(svector<int>& roots, unsigned x) {
unsigned old_x = x;
while (roots[x] >= 0) {
x = roots[x];
}
SASSERT(roots[x] < 0);
if (old_x != x) {
roots[old_x] = x;
}
return x;
}
static void merge(svector<int>& roots, unsigned x, unsigned y) {
x = find(roots, x);
y = find(roots, y);
SASSERT(roots[x] < 0 && roots[y] < 0);
if (x == y) {
return;
}
if (roots[x] > roots[y]) {
std::swap(x, y);
}
SASSERT(roots[x] <= roots[y]);
roots[x] += roots[y];
roots[y] = x;
}
template<typename Ext>
dl_var network_flow<Ext>::get_final(int start) {
int n = start;
while (m_depth[m_thread[n]] > m_depth[start]) {
n = m_thread[n];
}
return n;
}
template<typename Ext> template<typename Ext>
bool network_flow<Ext>::edge_in_tree(edge_id id) const { bool network_flow<Ext>::edge_in_tree(edge_id id) const {
return m_states[id] == BASIS; return m_states[id] == BASIS;
@ -540,202 +300,56 @@ namespace smt {
return edge_in_tree(get_edge_id(src,dst)); return edge_in_tree(get_edge_id(src,dst));
} }
/**
\brief Check invariants of main data-structures.
Spanning tree of m_graph + root is represented using:
svector<edge_state> m_states; edge_id |-> edge_state
svector<bool> m_upwards; node |-> bool
svector<node> m_pred; node |-> node
svector<int> m_depth; node |-> int
svector<node> m_thread; node |-> node
Tree is determined by m_pred:
- m_pred[root] == -1
- m_pred[n] = m != n for each node n, acyclic until reaching root.
- m_depth[m_pred[n]] + 1 == m_depth[n] for each n != root
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<typename Ext> template<typename Ext>
bool network_flow<Ext>::check_well_formed() { bool network_flow<Ext>::check_well_formed() {
node root = m_pred.size()-1; SASSERT(tree.check_well_formed());
// Check that m_thread traverses each node.
// This gets checked using union-find as well.
svector<bool> found(m_thread.size(), false);
found[root] = true;
for (node x = m_thread[root]; x != root; x = m_thread[x]) {
SASSERT(x != m_thread[x]);
found[x] = true;
}
for (unsigned i = 0; i < found.size(); ++i) {
SASSERT(found[i]);
}
// m_pred is acyclic, and points to root.
SASSERT(m_pred[root] == -1);
SASSERT(m_depth[root] == 0);
for (node i = 0; i < root; ++i) {
SASSERT(m_depth[m_pred[i]] < m_depth[i]);
}
// m_upwards show correct direction
for (unsigned i = 0; i < m_upwards.size(); ++i) {
node p = m_pred[i];
edge_id id;
SASSERT(!m_upwards[i] || m_graph.get_edge_id(i, p, id));
}
// m_depth[x] denotes distance from x to the root node
for (node x = m_thread[root]; x != root; x = m_thread[x]) {
SASSERT(m_depth[x] > 0);
SASSERT(m_depth[x] == m_depth[m_pred[x]] + 1);
}
// m_thread forms a spanning tree over [0..root]
// Union-find structure
svector<int> roots(m_pred.size(), -1);
for (node x = m_thread[root]; x != root; x = m_thread[x]) {
node y = m_pred[x];
// We are now going to check the edge between x and y
SASSERT(find(roots, x) != find(roots, y));
merge(roots, x, y);
}
// All nodes belong to the same spanning tree
for (unsigned i = 0; i < roots.size(); ++i) {
SASSERT(roots[i] + roots.size() == 0 || roots[i] >= 0);
}
// m_flows are zero on non-basic edges // m_flows are zero on non-basic edges
for (unsigned i = 0; i < m_flows.size(); ++i) { for (unsigned i = 0; i < m_flows.size(); ++i) {
SASSERT(m_states[i] == BASIS || m_flows[i].is_zero()); SASSERT(m_states[i] == BASIS || m_flows[i].is_zero());
} }
// m_upwards show correct direction
for (unsigned i = 0; i < m_potentials.size(); ++i) {
node p = tree.get_parent(i);
edge_id id;
SASSERT(!tree.get_arc_direction(i) || m_graph.get_edge_id(i, p, id));
}
return true; return true;
} }
template<typename Ext> template<typename Ext>
bool network_flow<Ext>::is_preorder_traversal(node start, node end) { std::string network_flow<Ext>:: display_spanning_tree() {
++m_step;;
// get children of start: std::ostringstream oss;
uint_set children; std::string prefix = "T";
children.insert(start); prefix.append(std::to_string(m_step));
node root = m_pred.size()-1; prefix.append("_");
for (int i = 0; i < root; ++i) { unsigned root = m_graph.get_num_nodes() - 1;
for (int j = 0; j < root; ++j) { for (unsigned i = 0; i < root; ++i) {
if (children.contains(m_pred[j])) { oss << prefix << i << "[shape=circle,label=\"" << prefix << i << " [";
children.insert(j); oss << m_potentials[i] << "/" << m_balances[i] << "]\"];\n";
}
}
}
// visit children using m_thread
children.remove(start);
do {
start = m_thread[start];
SASSERT(children.contains(start));
children.remove(start);
}
while (start != end);
SASSERT(children.empty());
return true;
} }
oss << prefix << root << "[shape=doublecircle,label=\"" << prefix << root << " [";
oss << m_potentials[root] << "/" << m_balances[root] << "]\"];\n";
} vector<edge> const & es = m_graph.get_all_edges();
for (unsigned i = 0; i < es.size(); ++i) {
#endif edge const & e = es[i];
oss << prefix << e.get_source() << " -> " << prefix << e.get_target();
#if 0 if (m_states[i] == BASIS) {
oss << "[color=red,penwidth=3.0,label=\"" << m_flows[i] << "/" << e.get_weight() << "\"];\n";
// At this point m_pred and m_upwards have been updated.
TRACE("network_flow", tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Upwards", m_upwards););
//
node x = m_final[p];
node y = m_thread[x];
node z = m_final[q];
// ----
// update m_final.
// Do this before updating data structures
node gamma_p = m_pred[m_thread[m_final[p]]];
node gamma_v = m_pred[m_thread[m_final[v]]];
node theta = m_thread[m_final[v]];
// Check that f(u) is not in T_v
bool found_final_u = false;
for (node n = v; n != theta; n = m_thread[n]) {
if (n == m_final[u]) {
found_final_u = true;
break;
}
}
node phi = find_rev_thread(v);
node delta = found_final_u ? phi : m_final[u];
TRACE("network_flow", tout << "Graft T_q and T_r'\n";);
for (node n = p; n != gamma_p && n != -1; n = m_pred[n]) {
TRACE("network_flow", tout << "1: m_final[" << n << "] |-> " << z << "\n";);
m_final[n] = z;
}
//
m_thread[x] = q;
m_thread[z] = y;
TRACE("network_flow", tout << "Update T_r'\n";);
m_thread[phi] = theta;
for (node n = u; n != gamma_v && n != -1; n = m_pred[n]) {
TRACE("network_flow", tout << "2: m_final[" << n << "] |-> " << delta << "\n";);
m_final[n] = delta;
}
TRACE("network_flow", tout << pp_vector("Last Successors", m_final, true) << pp_vector("Depths", m_depth););
if (v != q) {
TRACE("network_flow", tout << "Reroot T_v at q\n";);
node alpha1, alpha2;
node prev = q;
for (node n = v; n != q && n != -1; n = m_pred[n]) {
// Find all immediate successors of n
node t1 = m_thread[n];
node t2 = m_thread[m_final[t1]];
node t3 = m_thread[m_final[t2]];
if (t1 == m_pred[n]) {
alpha1 = t2;
alpha2 = t3;
}
else if (t2 == m_pred[n]) {
alpha1 = t1;
alpha2 = t3;
} }
else { else {
alpha1 = t1; oss << "[label=\"" << m_flows[i] << "/" << e.get_weight() << "\"];\n";
alpha2 = t2;
} }
m_thread[n] = alpha1;
m_thread[m_final[alpha1]] = alpha2;
m_thread[m_final[alpha2]] = prev;
prev = n;
} }
m_thread[m_final[q]] = prev; oss << std::endl;
return oss.str();
} }
}
#endif #endif

356
src/smt/spanning_tree.cpp Normal file
View file

@ -0,0 +1,356 @@
/*++
Copyright (c) 2013 Microsoft Corporation
Module Name:
spanning_tree.cpp
Abstract:
Author:
Anh-Dung Phan (t-anphan) 2013-11-06
Notes:
--*/
#include <sstream>
#include "spanning_tree.h"
#include "debug.h"
#include "vector.h"
#include "uint_set.h"
#include "trace.h"
namespace smt {
/**
swap v and q in tree.
- fixup m_thread
- fixup m_pred
Case 1: final(q) == final(v)
-------
Old thread: prev -> v -*-> alpha -> q -*-> final(q) -> next
New thread: prev -> q -*-> final(q) -> v -*-> alpha -> next
Case 2: final(q) != final(v)
-------
Old thread: prev -> v -*-> alpha -> q -*-> final(q) -> beta -*-> final(v) -> next
New thread: prev -> q -*-> final(q) -> v -*-> alpha -> beta -*-> final(v) -> next
*/
void thread_spanning_tree::swap_order(node q, node v) {
SASSERT(q != v);
SASSERT(m_pred[q] == v);
SASSERT(is_preorder_traversal(v, get_final(v)));
node prev = find_rev_thread(v);
node final_q = get_final(q);
node final_v = get_final(v);
node next = m_thread[final_v];
node alpha = find_rev_thread(q);
if (final_q == final_v) {
m_thread[final_q] = v;
m_thread[alpha] = next;
}
else {
node beta = m_thread[final_q];
m_thread[final_q] = v;
m_thread[alpha] = beta;
}
m_thread[prev] = q;
m_pred[v] = q;
SASSERT(is_preorder_traversal(q, get_final(q)));
}
/**
\brief find node that points to 'n' in m_thread
*/
node thread_spanning_tree::find_rev_thread(node n) const {
node ancestor = m_pred[n];
SASSERT(ancestor != -1);
while (m_thread[ancestor] != n) {
ancestor = m_thread[ancestor];
}
return ancestor;
}
void thread_spanning_tree::fix_depth(node start, node end) {
SASSERT(m_pred[start] != -1);
m_depth[start] = m_depth[m_pred[start]]+1;
while (start != end) {
start = m_thread[start];
m_depth[start] = m_depth[m_pred[start]]+1;
}
}
node thread_spanning_tree::get_final(int start) {
int n = start;
while (m_depth[m_thread[n]] > m_depth[start]) {
n = m_thread[n];
}
return n;
}
bool thread_spanning_tree::is_preorder_traversal(node start, node end) {
// get children of start
uint_set children;
children.insert(start);
node root = m_pred.size()-1;
for (int i = 0; i < root; ++i) {
for (int j = 0; j < root; ++j) {
if (children.contains(m_pred[j])) {
children.insert(j);
}
}
}
// visit children using m_thread
children.remove(start);
do {
start = m_thread[start];
SASSERT(children.contains(start));
children.remove(start);
}
while (start != end);
SASSERT(children.empty());
return true;
}
bool thread_spanning_tree::is_ancestor_of(node ancestor, node child) {
for (node n = child; n != -1; n = m_pred[n]) {
if (n == ancestor) {
return true;
}
}
return false;
}
static unsigned find(svector<int>& roots, unsigned x) {
unsigned old_x = x;
while (roots[x] >= 0) {
x = roots[x];
}
SASSERT(roots[x] < 0);
if (old_x != x) {
roots[old_x] = x;
}
return x;
}
static void merge(svector<int>& roots, unsigned x, unsigned y) {
x = find(roots, x);
y = find(roots, y);
SASSERT(roots[x] < 0 && roots[y] < 0);
if (x == y) {
return;
}
if (roots[x] > roots[y]) {
std::swap(x, y);
}
SASSERT(roots[x] <= roots[y]);
roots[x] += roots[y];
roots[y] = x;
}
void thread_spanning_tree::initialize(svector<bool> const & upwards, int num_nodes) {
m_pred.resize(num_nodes + 1);
m_depth.resize(num_nodes + 1);
m_thread.resize(num_nodes + 1);
m_upwards.resize(num_nodes + 1);
node root = num_nodes;
m_pred[root] = -1;
m_depth[root] = 0;
m_thread[root] = 0;
// Create artificial edges from/to root node to/from other nodes and initialize the spanning tree
for (int i = 0; i < num_nodes; ++i) {
m_pred[i] = root;
m_depth[i] = 1;
m_thread[i] = i + 1;
m_upwards[i] = upwards[i];
}
TRACE("network_flow", {
tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread);
tout << pp_vector("Depths", m_depth) << pp_vector("Upwards", m_upwards);
});
}
node thread_spanning_tree::get_common_ancestor(node u, node v) {
while (u != v) {
if (m_depth[u] > m_depth[v])
u = m_pred[u];
else
v = m_pred[v];
}
return u;
}
void thread_spanning_tree::get_descendants(node start, svector<node>& descendants) {
descendants.reset();
node u = start;
while (m_depth[m_thread[u]] > m_depth[start]) {
descendants.push_back(u);
u = m_thread[u];
}
}
void thread_spanning_tree::get_ancestors(node start, svector<node>& ancestors) {
ancestors.reset();
while (m_pred[start] != -1) {
ancestors.push_back(start);
start = m_pred[start];
}
}
/**
\brief add entering_edge, remove leaving_edge from spanning tree.
Old tree: New tree:
root root
/ \ / \
x y x y
/ \ / \ / \ / \
u s u s
| / /
v w v w
/ \ \ / \ \
z p z p
\ \ /
q q
*/
void thread_spanning_tree::update(node p, node q, node u, node v) {
bool q_upwards = false;
// v is parent of u so T_u does not contain root node
if (m_pred[u] == v) {
std::swap(u, v);
}
SASSERT(m_pred[v] == u);
if (is_ancestor_of(v, p)) {
std::swap(p, q);
q_upwards = true;
}
SASSERT(is_ancestor_of(v, q));
TRACE("network_flow", {
tout << "update_spanning_tree: (" << p << ", " << q << ") enters, (";
tout << u << ", " << v << ") leaves\n";
});
// Update m_pred (for nodes in the stem from q to v)
// Note: m_pred[v] == u
// Initialize m_upwards[q] = q_upwards
bool prev_upwards = q_upwards;
node old_pred = m_pred[q];
if (q != v) {
for (node n = q; n != u; ) {
SASSERT(old_pred != u || n == v); // the last processed node is v
TRACE("network_flow", {
tout << pp_vector("Predecessors", m_pred, true);
});
SASSERT(-1 != m_pred[old_pred]);
int next_old_pred = m_pred[old_pred];
swap_order(n, old_pred);
std::swap(m_upwards[n], prev_upwards);
prev_upwards = !prev_upwards; // flip previous version of upwards.
n = old_pred;
old_pred = next_old_pred;
}
}
m_pred[q] = p;
// m_thread were updated.
// update the depth.
fix_depth(q, get_final(q));
TRACE("network_flow", {
tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread);
tout << pp_vector("Depths", m_depth) << pp_vector("Upwards", m_upwards);
});
}
/**
\brief Check invariants of main data-structures.
Spanning tree of m_graph + root is represented using:
svector<edge_state> m_states; edge_id |-> edge_state
svector<bool> m_upwards; node |-> bool
svector<node> m_pred; node |-> node
svector<int> m_depth; node |-> int
svector<node> m_thread; node |-> node
Tree is determined by m_pred:
- m_pred[root] == -1
- m_pred[n] = m != n for each node n, acyclic until reaching root.
- m_depth[m_pred[n]] + 1 == m_depth[n] for each n != root
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
*/
bool thread_spanning_tree::check_well_formed() {
node root = m_pred.size()-1;
// Check that m_thread traverses each node.
// This gets checked using union-find as well.
svector<bool> found(m_thread.size(), false);
found[root] = true;
for (node x = m_thread[root]; x != root; x = m_thread[x]) {
SASSERT(x != m_thread[x]);
found[x] = true;
}
for (unsigned i = 0; i < found.size(); ++i) {
SASSERT(found[i]);
}
// m_pred is acyclic, and points to root.
SASSERT(m_pred[root] == -1);
SASSERT(m_depth[root] == 0);
for (node i = 0; i < root; ++i) {
SASSERT(m_depth[m_pred[i]] < m_depth[i]);
}
// m_depth[x] denotes distance from x to the root node
for (node x = m_thread[root]; x != root; x = m_thread[x]) {
SASSERT(m_depth[x] > 0);
SASSERT(m_depth[x] == m_depth[m_pred[x]] + 1);
}
// m_thread forms a spanning tree over [0..root]
// Union-find structure
svector<int> roots(m_pred.size(), -1);
for (node x = m_thread[root]; x != root; x = m_thread[x]) {
node y = m_pred[x];
// We are now going to check the edge between x and y
SASSERT(find(roots, x) != find(roots, y));
merge(roots, x, y);
}
// All nodes belong to the same spanning tree
for (unsigned i = 0; i < roots.size(); ++i) {
SASSERT(roots[i] + roots.size() == 0 || roots[i] >= 0);
}
return true;
}
bool thread_spanning_tree::get_arc_direction(node start) const {
return m_upwards[start];
}
node thread_spanning_tree::get_parent(node start) {
return m_pred[start];
}
}

63
src/smt/spanning_tree.h Normal file
View file

@ -0,0 +1,63 @@
/*++
Copyright (c) 2013 Microsoft Corporation
Module Name:
spanning_tree.h
Abstract:
Represent spanning trees with needed operations for Network Simplex
Author:
Anh-Dung Phan (t-anphan) 2013-11-06
Notes:
--*/
#ifndef _SPANNING_TREE_H_
#define _SPANNING_TREE_H_
#include "spanning_tree_base.h"
namespace smt {
class thread_spanning_tree : virtual public spanning_tree_base {
private:
// Store the parent of a node i in the spanning tree
svector<node> m_pred;
// Store the number of edge on the path from node i to the root
svector<int> m_depth;
// Store the pointer from node i to the next node in depth-first search order
svector<node> m_thread;
// m_upwards[i] is true if the corresponding edge
// (i, m_pred[i]) points upwards (pointing toward the root node)
svector<bool> m_upwards;
void swap_order(node q, node v);
node find_rev_thread(node n) const;
void fix_depth(node start, node end);
node get_final(int start);
bool is_preorder_traversal(node start, node end);
bool is_ancestor_of(node ancestor, node child);
public:
void initialize(svector<bool> const & upwards, int num_nodes);
void get_descendants(node start, svector<node>& descendants);
void get_ancestors(node start, svector<node>& ancestors);
node get_common_ancestor(node u, node v);
void update(node p, node q, node u, node v);
bool check_well_formed();
// TODO: remove these two unnatural functions
bool get_arc_direction(node start) const;
node get_parent(node start);
};
}
#endif

71
src/smt/spanning_tree_base.h Normal file
View file

@ -0,0 +1,71 @@
/*++
Copyright (c) 2013 Microsoft Corporation
Module Name:
spanning_tree_base.h
Abstract:
Represent spanning trees with needed operations for Network Simplex
Author:
Anh-Dung Phan (t-anphan) 2013-11-06
Notes:
--*/
#ifndef _SPANNING_TREE_BASE_H_
#define _SPANNING_TREE_BASE_H_
#include "util.h"
#include "vector.h"
namespace smt {
typedef int node;
template<typename TV>
inline std::string pp_vector(std::string const & label, TV v, bool has_header = false) {
std::ostringstream oss;
if (has_header) {
oss << "Index ";
for (unsigned i = 0; i < v.size(); ++i) {
oss << i << " ";
}
oss << std::endl;
}
oss << label << " ";
for (unsigned i = 0; i < v.size(); ++i) {
oss << v[i] << " ";
}
oss << std::endl;
return oss.str();
}
class spanning_tree_base {
public:
spanning_tree_base() {};
virtual ~spanning_tree_base() {};
virtual void initialize(svector<bool> const & upwards, int num_nodes) = 0;
/**
\brief Get all descendants of a node including itself
*/
virtual void get_descendants(node start, svector<node>& descendants) = 0;
/**
\brief Get all ancestors of a node including itself
*/
virtual void get_ancestors(node start, svector<node>& ancestors) = 0;
virtual node get_common_ancestor(node u, node v) = 0;
virtual void update(node p, node q, node u, node v) = 0;
virtual bool check_well_formed() = 0;
// TODO: remove these two unnatural functions
virtual bool get_arc_direction(node start) const = 0;
virtual node get_parent(node start) = 0;
};
}
#endif