diff --git a/src/smt/network_flow.h b/src/smt/network_flow.h index a90e901e8..489cf351e 100644 --- a/src/smt/network_flow.h +++ b/src/smt/network_flow.h @@ -32,12 +32,10 @@ Notes: #include"inf_rational.h" #include"diff_logic.h" +#include"spanning_tree.h" namespace smt { - template - std::string pp_vector(std::string const & label, TV v, bool has_header = false); - // Solve minimum cost flow problem using Network Simplex algorithm template class network_flow : private Ext { @@ -51,7 +49,9 @@ namespace smt { typedef dl_graph graph; typedef typename Ext::numeral numeral; typedef typename Ext::fin_numeral fin_numeral; + graph m_graph; + thread_spanning_tree tree; // Denote supply/demand b_i on node i vector m_balances; @@ -67,18 +67,7 @@ namespace smt { svector m_states; - // m_upwards[i] is true if the corresponding edge - // (i, m_pred[i]) points upwards (pointing toward the root node) - svector m_upwards; - - // Store the parent of a node i in the spanning tree - svector m_pred; - // Store the number of edge on the path from node i to the root - svector m_depth; - // Store the pointer from node i to the next node in depth-first search order - svector m_thread; - // Store a final node of the sub tree rooted at node i - svector m_final; + unsigned m_step; edge_id m_entering_edge; edge_id m_leaving_edge; @@ -86,14 +75,11 @@ namespace smt { optional m_delta; bool m_in_edge_dir; - unsigned m_step; - // Initialize the network with a feasible spanning tree void initialize(); edge_id get_edge_id(dl_var source, dl_var target) const; - void update_potentials(); void update_flows(); @@ -113,22 +99,8 @@ namespace smt { bool edge_in_tree(edge_id id) 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 is_preorder_traversal(node start, node end); - node get_final(int start); - public: network_flow(graph & g, vector const & balances); diff --git a/src/smt/network_flow_def.h b/src/smt/network_flow_def.h index 96557301f..834cfc075 100644 --- a/src/smt/network_flow_def.h +++ b/src/smt/network_flow_def.h @@ -25,24 +25,6 @@ Notes: namespace smt { - template - 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 network_flow::network_flow(graph & g, vector const & balances) : m_balances(balances) { @@ -63,14 +45,9 @@ namespace smt { unsigned num_edges = m_graph.get_num_edges(); m_balances.resize(num_nodes); - m_potentials.resize(num_nodes); + m_potentials.resize(num_nodes); - m_upwards.resize(num_nodes); - m_pred.resize(num_nodes); - m_depth.resize(num_nodes); - m_thread.resize(num_nodes); - - m_step = 0; + tree = thread_spanning_tree(); } template @@ -79,13 +56,10 @@ namespace smt { // Create an artificial root node to construct initial spanning tree unsigned num_nodes = m_graph.get_num_nodes(); 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(); - m_graph.init_var(root); + m_potentials[root] = numeral::zero(); fin_numeral sum_supply = fin_numeral::zero(); for (unsigned i = 0; i < m_balances.size(); ++i) { @@ -99,31 +73,31 @@ namespace smt { m_states.fill(LOWER); // Create artificial edges from/to root node to/from other nodes and initialize the spanning tree + svector upwards(num_nodes, false); for (unsigned i = 0; i < num_nodes; ++i) { - m_upwards[i] = !m_balances[i].is_neg(); - m_pred[i] = root; - m_depth[i] = 1; - m_thread[i] = i + 1; + upwards[i] = !m_balances[i].is_neg(); m_states[num_edges + i] = BASIS; - node src = m_upwards[i] ? i : root; - node tgt = m_upwards[i] ? root : i; - m_flows[num_edges + i] = m_upwards[i] ? m_balances[i] : -m_balances[i]; + node src = upwards[i] ? i : root; + node tgt = upwards[i] ? root : i; + m_flows[num_edges + i] = upwards[i] ? m_balances[i] : -m_balances[i]; m_graph.add_edge(src, tgt, numeral::one(), explanation()); } + tree.initialize(upwards, num_nodes); + // Compute initial potentials - node u = m_thread[root]; - while (u != root) { - node v = m_pred[u]; + svector descendants; + tree.get_descendants(root, descendants); + // 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); - m_potentials[u] = m_potentials[v] + (m_upwards[u] ? - m_graph.get_weight(e_id) : m_graph.get_weight(e_id)); - u = m_thread[u]; + m_potentials[u] = m_potentials[v] + (tree.get_arc_direction(u) ? - m_graph.get_weight(e_id) : m_graph.get_weight(e_id)); } 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); - tout << pp_vector("Potentials", m_potentials) << pp_vector("Flows", m_flows); + tout << pp_vector("Potentials", m_potentials, true) << pp_vector("Flows", m_flows); }); TRACE("network_flow", tout << "Spanning tree:\n" << display_spanning_tree();); SASSERT(check_well_formed()); @@ -131,9 +105,9 @@ namespace smt { template edge_id network_flow::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; - 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; } @@ -143,9 +117,11 @@ namespace smt { node src = m_graph.get_source(m_entering_edge); node tgt = m_graph.get_target(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]); - node last = m_thread[get_final(src)]; - for (node u = src; u != last; u = m_thread[u]) { + numeral change = tree.get_arc_direction(src) ? (-cost - m_potentials[src] + m_potentials[tgt]) : (cost + m_potentials[src] - m_potentials[tgt]); + svector descendants; + tree.get_descendants(src, descendants); + for (unsigned i = 0; i < descendants.size(); ++i) { + node u = descendants[i]; m_potentials[u] += change; } 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); m_flows[m_entering_edge] += val; node source = m_graph.get_source(m_entering_edge); - for (unsigned u = source; u != m_join_node; u = m_pred[u]) { - edge_id e_id = get_edge_id(u, m_pred[u]); - m_flows[e_id] += m_upwards[u] ? -val : val; + svector ancestors; + tree.get_ancestors(source, ancestors); + 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); - for (unsigned u = target; u != m_join_node; u = m_pred[u]) { - edge_id e_id = get_edge_id(u, m_pred[u]); - m_flows[e_id] += m_upwards[u] ? val : -val; + tree.get_ancestors(target,ancestors); + 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; } TRACE("network_flow", tout << pp_vector("Flows", m_flows, true);); } @@ -202,41 +184,37 @@ namespace smt { if (m_states[m_entering_edge] == UPPER) { std::swap(source, target); } - node u = source, v = target; - while (u != v) { - 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; + + m_join_node = tree.get_common_ancestor(source, target); + TRACE("network_flow", tout << "Found join node " << m_join_node << std::endl;); m_delta.set_invalid(); node src, tgt; + // Send flows along the path from source to the ancestor - for (unsigned u = source; u != m_join_node; u = m_pred[u]) { - edge_id e_id = get_edge_id(u, m_pred[u]); - if (m_upwards[u] && (!m_delta || m_flows[e_id] < *m_delta)) { + svector ancestors; + tree.get_ancestors(source, ancestors); + 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]; src = u; - tgt = m_pred[u]; + tgt = tree.get_parent(u); SASSERT(edge_in_tree(src,tgt)); m_in_edge_dir = true; } } // Send flows along the path from target to the ancestor - for (unsigned u = target; u != m_join_node; u = m_pred[u]) { - edge_id e_id = get_edge_id(u, m_pred[u]); - if (!m_upwards[u] && (!m_delta || m_flows[e_id] <= *m_delta)) { + tree.get_ancestors(target, ancestors); + 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]; src = u; - tgt = m_pred[u]; + tgt = tree.get_parent(u); SASSERT(edge_in_tree(src,tgt)); m_in_edge_dir = false; } @@ -254,196 +232,14 @@ namespace smt { return false; } - template - bool network_flow::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 - dl_var network_flow::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 - void network_flow::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 void network_flow::update_spanning_tree() { node p = m_graph.get_source(m_entering_edge); node q = m_graph.get_target(m_entering_edge); node u = m_graph.get_source(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 - 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 - void network_flow::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 - std::string network_flow::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 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(); + + tree.update(p, q, u, v); } // Minimize cost flows @@ -494,42 +290,6 @@ namespace smt { return m_objective_value; } - static unsigned find(svector& 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& 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 - dl_var network_flow::get_final(int start) { - int n = start; - while (m_depth[m_thread[n]] > m_depth[start]) { - n = m_thread[n]; - } - return n; - } - template bool network_flow::edge_in_tree(edge_id id) const { return m_states[id] == BASIS; @@ -540,202 +300,56 @@ namespace smt { 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 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 - - 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 bool network_flow::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 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 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); - } + SASSERT(tree.check_well_formed()); // m_flows are zero on non-basic edges for (unsigned i = 0; i < m_flows.size(); ++i) { 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; } template - bool network_flow::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); + std::string network_flow:: 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 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"; } } - } - // 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 << std::endl; + return oss.str(); } - } #endif - -#if 0 - - // 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 { - alpha1 = t1; - 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; - } -#endif diff --git a/src/smt/spanning_tree.cpp b/src/smt/spanning_tree.cpp new file mode 100644 index 000000000..21c3fcfb4 --- /dev/null +++ b/src/smt/spanning_tree.cpp @@ -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 +#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& 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& 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 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& 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& 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 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 + + 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 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 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]; + } +} \ No newline at end of file diff --git a/src/smt/spanning_tree.h b/src/smt/spanning_tree.h new file mode 100644 index 000000000..6769a6a18 --- /dev/null +++ b/src/smt/spanning_tree.h @@ -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 m_pred; + // Store the number of edge on the path from node i to the root + svector m_depth; + // Store the pointer from node i to the next node in depth-first search order + svector m_thread; + + // m_upwards[i] is true if the corresponding edge + // (i, m_pred[i]) points upwards (pointing toward the root node) + svector 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 const & upwards, int num_nodes); + void get_descendants(node start, svector& descendants); + void get_ancestors(node start, svector& 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 diff --git a/src/smt/spanning_tree_base.h b/src/smt/spanning_tree_base.h new file mode 100644 index 000000000..6be090711 --- /dev/null +++ b/src/smt/spanning_tree_base.h @@ -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 + 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 const & upwards, int num_nodes) = 0; + /** + \brief Get all descendants of a node including itself + */ + virtual void get_descendants(node start, svector& descendants) = 0; + /** + \brief Get all ancestors of a node including itself + */ + virtual void get_ancestors(node start, svector& 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 \ No newline at end of file