From 905f230b8fb2b6d249c324344efe6aac82924b6e Mon Sep 17 00:00:00 2001 From: Anh-Dung Phan Date: Tue, 29 Oct 2013 14:20:29 -0700 Subject: [PATCH 1/2] Add pretty printing for network_flow Reuse the original graph as much as possible --- src/smt/diff_logic.h | 2 + src/smt/network_flow.h | 27 ++--- src/smt/network_flow_def.h | 203 ++++++++++++++++++++++++++------ src/smt/theory_diff_logic_def.h | 25 ++-- 4 files changed, 190 insertions(+), 67 deletions(-) diff --git a/src/smt/diff_logic.h b/src/smt/diff_logic.h index 71b638b9e..b96a038eb 100644 --- a/src/smt/diff_logic.h +++ b/src/smt/diff_logic.h @@ -281,6 +281,8 @@ public: unsigned get_num_edges() const { return m_edges.size(); } + unsigned get_num_nodes() const { return m_out_edges.size(); } + dl_var get_source(edge_id id) const { return m_edges[id].get_source(); } dl_var get_target(edge_id id) const { return m_edges[id].get_target(); } diff --git a/src/smt/network_flow.h b/src/smt/network_flow.h index 37afda136..b0539811e 100644 --- a/src/smt/network_flow.h +++ b/src/smt/network_flow.h @@ -47,7 +47,7 @@ namespace smt { typedef dl_graph graph; typedef typename Ext::numeral numeral; typedef typename Ext::fin_numeral fin_numeral; - graph m_graph; + graph & m_graph; // Denote supply/demand b_i on node i vector m_balances; @@ -57,10 +57,7 @@ namespace smt { // Keep optimal solution of the min cost flow problem numeral m_objective_value; - - // Costs on edges - vector m_costs; - + // Basic feasible flows vector m_flows; @@ -79,18 +76,12 @@ namespace smt { svector m_rev_thread; // Store a final node of the sub tree rooted at node i svector m_final; - // Number of nodes in the sub tree rooted at node i - svector m_num_node; edge_id m_entering_edge; edge_id m_leaving_edge; node m_join_node; numeral m_delta; - public: - - network_flow(graph & g, vector const & balances); - // Initialize the network with a feasible spanning tree void initialize(); @@ -107,13 +98,21 @@ namespace smt { bool choose_leaving_edge(); void update_spanning_tree(); - - // Compute the optimal solution - numeral get_optimal_solution(vector & result, bool is_dual); + std::string pp_vector(std::string const & label, svector v, bool has_header = false); + std::string pp_vector(std::string const & label, vector v, bool has_header = false); + + public: + + network_flow(graph & g, vector const & balances); + // Minimize cost flows // Return true if found an optimal solution, and return false if unbounded bool min_cost(); + + // 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 f39042bac..07548d2f6 100644 --- a/src/smt/network_flow_def.h +++ b/src/smt/network_flow_def.h @@ -25,25 +25,15 @@ Notes: namespace smt { template - network_flow::network_flow(graph & g, vector const& balances) : + network_flow::network_flow(graph & g, vector const & balances) : m_graph(g), m_balances(balances) { - unsigned num_nodes = m_balances.size() + 1; + unsigned num_nodes = m_graph.get_num_nodes() + 1; unsigned num_edges = m_graph.get_num_edges(); - vector const & es = m_graph.get_all_edges(); - for (unsigned i = 0; i < num_edges; ++i) { - fin_numeral cost(es[i].get_weight().get_rational()); - m_costs.push_back(cost); - } - m_balances.resize(num_nodes); - for (unsigned i = 0; i < balances.size(); ++i) { - m_costs.push_back(balances[i]); - } - m_potentials.resize(num_nodes); - m_costs.resize(num_edges); + m_potentials.resize(num_nodes); m_flows.resize(num_edges); m_states.resize(num_edges); @@ -53,19 +43,18 @@ namespace smt { m_thread.resize(num_nodes); m_rev_thread.resize(num_nodes); m_final.resize(num_nodes); - m_num_node.resize(num_nodes); } template void network_flow::initialize() { + TRACE("network_flow", tout << "initialize...\n";); // Create an artificial root node to construct initial spanning tree - unsigned num_nodes = m_balances.size(); + 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_thread[root] = 0; m_rev_thread[0] = root; - m_num_node[root] = num_nodes + 1; m_final[root] = root - 1; m_potentials[root] = numeral::zero(); @@ -85,14 +74,17 @@ namespace smt { m_depth[i] = 1; m_thread[i] = i + 1; m_final[i] = i; - m_rev_thread[i] = (i = 0) ? root : i - 1; - m_num_node[i] = 1; + m_rev_thread[i] = (i == 0) ? root : i - 1; m_states[num_edges + i] = BASIS; } + + TRACE("network_flow", tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread) + << pp_vector("Reverse Threads", m_rev_thread) << pp_vector("Last Successors", m_final) << pp_vector("Depths", m_depth);); } template void network_flow::update_potentials() { + TRACE("network_flow", tout << "update_potentials...\n";); node src = m_graph.get_source(m_entering_edge); node tgt = m_graph.get_source(m_entering_edge); numeral cost = m_graph.get_weight(m_entering_edge); @@ -102,31 +94,37 @@ namespace smt { for (node u = src; u != last; u = m_thread[u]) { m_potentials[u] += change; } + TRACE("network_flow", tout << pp_vector("Potentials", m_potentials, true);); } template void network_flow::update_flows() { + TRACE("network_flow", tout << "update_flows...\n";); numeral val = m_state[m_entering_edge] == NON_BASIS ? numeral::zero() : m_delta; m_flows[m_entering_edge] += val; - for (unsigned u = m_graph.get_source(m_entering_edge); u != m_join_node; u = m_pred[u]) { + 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; m_graph.get_edge_id(u, m_pred[u], e_id); m_flows[e_id] += m_upwards[u] ? -val : val; } - for (unsigned u = m_graph.get_target(m_entering_edge); u != m_join_node; u = m_pred[u]) { + 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; m_graph.get_edge_id(u, m_pred[u], e_id); m_flows[e_id] += m_upwards[u] ? val : -val; } + 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";); vector const & es = m_graph.get_all_edges(); for (unsigned int i = 0; i < es.size(); ++i) { edge const & e = es[i]; edge_id e_id; - if (e.is_enabled() && m_graph.get_edge_id(e.get_source(), e.get_target(), e_id) && m_states[e_id] == BASIS) { + if (e.is_enabled() && m_graph.get_edge_id(e.get_source(), e.get_target(), e_id) && m_states[e_id] == NON_BASIS) { node source = e.get_source(); node target = e.get_target(); numeral cost = e.get_weight() - m_potentials[source] + m_potentials[target]; @@ -134,15 +132,18 @@ namespace smt { // TODO: add multiple pivoting strategies if (cost < numeral::zero()) { m_entering_edge = e_id; + TRACE("network_flow", tout << "Found entering edge " << e_id << " between node " << source << " and node " << target << "...\n";); return true; } } } + TRACE("network_flow", tout << "Found no entering edge... It's probably optimal.\n";); return false; } template bool network_flow::choose_leaving_edge() { + TRACE("network_flow", tout << "choose_leaving_edge...\n";); node source = m_graph.get_source(m_entering_edge); node target = m_graph.get_target(m_entering_edge); node u = source, v = target; @@ -158,9 +159,9 @@ namespace smt { } // Found first common ancestor of source and target m_join_node = u; - // FIXME: need to get truly finite value + // FIXME: need to get truly infinite value numeral infty = numeral(INT_MAX); - m_delta = infty; + m_delta = infty; 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]) { @@ -188,31 +189,143 @@ namespace smt { if (m_delta < infty) { m_graph.get_edge_id(src, tgt, m_leaving_edge); + TRACE("network_flow", tout << "Found leaving edge" << m_leaving_edge << "between node " << src << " and node " << tgt << "...\n";); return true; } + TRACE("network_flow", tout << "Can't find a leaving edge... The problem is unbounded.\n";); return false; } template - void network_flow::update_spanning_tree() { + void network_flow::update_spanning_tree() { + node src_in = m_graph.get_source(m_entering_edge); + node tgt_in = m_graph.get_target(m_entering_edge); + node src_out = m_graph.get_source(m_leaving_edge); + node tgt_out = m_graph.get_target(m_leaving_edge); + TRACE("network_flow", tout << "update_spanning_tree: (" << src_in << ", " << tgt_in << ") enters, (" + << src_out << ", " << tgt_out << ") leaves\n";); + node root = m_graph.get_num_nodes(); + node rev_thread_out = m_rev_thread[src_out]; + + node x = m_final[src_in]; + node y = m_thread[x]; + node z = m_final[src_out]; + + // Update m_pred (for nodes in the stem from tgt_in to tgt_out) + node u = tgt_in; + node last = m_pred[tgt_out]; + node parent = src_in; + while (u != last) { + node next = m_pred[u]; + m_pred[u] = parent; + u = next; + } + + // Graft T_q and T_r' + m_thread[x] = src_out; + m_thread[z] = y; + u = src_out; + while (u != m_final[src_out]) { + m_depth[u] += 1 + m_depth[src_in]; + u = m_pred[u]; + } + + node gamma = m_thread[m_final[src_in]]; + last = m_pred[gamma] != 0 ? gamma : root; + for (node u = src_in; u == last; u = m_pred[u]) { + m_final[u] = z; + } + + // Update T_r' + node phi = m_thread[tgt_out]; + node theta = m_thread[m_final[tgt_out]]; + m_thread[phi] = theta; + + gamma = m_thread[m_final[tgt_out]]; + // REVIEW: check f(u) is not in T_v + node delta = m_final[src_out] != m_final[tgt_out] ? m_final[src_out] : m_rev_thread[tgt_out]; + last = m_pred[gamma] != 0 ? gamma : root; + for (node u = src_in; u == last; u = m_pred[u]) { + m_final[u] = delta; + } + + // Reroot T_v at q + if (src_out != tgt_in) { + node u = m_pred[src_out]; + m_thread[m_final[src_out]] = u; + last = tgt_in; + node alpha1, alpha2; + unsigned count = 0; + while (u != last) { + // Find all immediate successors of u + node t1 = m_thread[u]; + node t2 = m_thread[m_final[t1]]; + node t3 = m_thread[m_final[t2]]; + if (t1 = m_pred[u]) { + alpha1 = t2; + alpha2 = t3; + } + else if (t2 = m_pred[u]) { + alpha1 = t1; + alpha2 = t3; + } + else { + alpha1 = t1; + alpha2 = t2; + } + m_thread[u] = alpha1; + m_thread[m_final[alpha1]] = alpha2; + u = m_pred[u]; + m_thread[m_final[alpha2]] = u; + // Decrease depth of all children in the subtree + ++count; + int d = m_depth[u] - count; + for (node v = m_thread[u]; v == m_final[u]; v = m_thread[v]) { + m_depth[v] -= d; + } + } + m_thread[m_final[alpha2]] = src_out; + } + + TRACE("network_flow", tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread) + << pp_vector("Reverse Threads", m_rev_thread) << pp_vector("Last Successors", m_final) << pp_vector("Depths", m_depth);); + } + + template + std::string network_flow::pp_vector(std::string const & label, svector 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(); } - // Get the optimal solution template - typename network_flow::numeral network_flow::get_optimal_solution(vector & result, bool is_dual) { - m_objective_value = numeral::zero(); - for (unsigned i = 0; i < m_flows.size(); ++i) { - m_objective_value += m_costs[i] * m_flows[i]; + std::string network_flow::pp_vector(std::string const & label, vector 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; } - result.reset(); - if (is_dual) { - result.append(m_potentials); + oss << label << " "; + for (unsigned i = 0; i < v.size(); ++i) { + oss << v[i] << " "; } - else { - result.append(m_flows); - } - return m_objective_value; + oss << std::endl; + return oss.str(); } // Minimize cost flows @@ -232,6 +345,24 @@ namespace smt { } return true; } + + // Get the optimal solution + template + typename network_flow::numeral network_flow::get_optimal_solution(vector & result, bool is_dual) { + m_objective_value = numeral::zero(); + for (unsigned i = 0; i < m_flows.size(); ++i) { + fin_numeral cost = m_graph.get_weight(i).get_rational(); + m_objective_value += cost * m_flows[i]; + } + result.reset(); + if (is_dual) { + result.append(m_potentials); + } + else { + result.append(m_flows); + } + return m_objective_value; + } } #endif diff --git a/src/smt/theory_diff_logic_def.h b/src/smt/theory_diff_logic_def.h index 9fa97097f..70fc681a7 100644 --- a/src/smt/theory_diff_logic_def.h +++ b/src/smt/theory_diff_logic_def.h @@ -1008,34 +1008,25 @@ bool theory_diff_logic::maximize(theory_var v) { verbose_stream() << "coefficient " << objective[i].second << " of theory_var " << objective[i].first << "\n"; } verbose_stream() << m_objective_consts[v] << "\n";); - NOT_IMPLEMENTED_YET(); - // Double the number of edges in the new graph - - // NSB review: what about disabled edges? They should not be added, right? - // For efficiency we probably want to reuse m_graph and keep extra edges on the side or add them to - // m_graph as well. - dl_graph g; - vector > const& es = m_graph.get_all_edges(); - for (unsigned i = 0; i < es.size(); ++i) { - dl_edge const & e = es[i]; - if (e.is_enabled()) { - g.enable_edge(g.add_edge(e.get_source(), e.get_target(), e.get_weight(), e.get_explanation())); - } - } // Objective coefficients now become balances - vector balances; + vector balances(m_graph.get_num_nodes()); + balances.fill(fin_numeral::zero()); for (unsigned i = 0; i < objective.size(); ++i) { fin_numeral balance(objective[i].second); - balances.push_back(balance); + balances[objective[i].first] = balance; } - network_flow net_flow(g, balances); + network_flow net_flow(m_graph, balances); bool is_optimal = net_flow.min_cost(); if (is_optimal) { vector potentials; m_objective_value = net_flow.get_optimal_solution(potentials, true); + std::cout << "Objective value of " << v << ": " << m_objective_value << std::endl; // TODO: return the model of the optimal solution from potential + } + else { + std::cout << "Unbounded" << std::endl; } return is_optimal; } From e715ccbc987b8fdafa9c80ecda95d2fcb3cbd015 Mon Sep 17 00:00:00 2001 From: Anh-Dung Phan Date: Tue, 29 Oct 2013 15:49:53 -0700 Subject: [PATCH 2/2] Minor updates --- src/smt/network_flow.h | 5 +++- src/smt/network_flow_def.h | 43 +++++++++++++++++++++------------ src/smt/theory_diff_logic_def.h | 2 +- 3 files changed, 33 insertions(+), 17 deletions(-) diff --git a/src/smt/network_flow.h b/src/smt/network_flow.h index b0539811e..5177af2e4 100644 --- a/src/smt/network_flow.h +++ b/src/smt/network_flow.h @@ -47,7 +47,7 @@ namespace smt { typedef dl_graph graph; typedef typename Ext::numeral numeral; typedef typename Ext::fin_numeral fin_numeral; - graph & m_graph; + graph m_graph; // Denote supply/demand b_i on node i vector m_balances; @@ -85,6 +85,9 @@ namespace smt { // Initialize the network with a feasible spanning tree void initialize(); + bool get_edge_id(dl_var source, dl_var target, edge_id & id); + + void update_potentials(); void update_flows(); diff --git a/src/smt/network_flow_def.h b/src/smt/network_flow_def.h index 07548d2f6..ba156c712 100644 --- a/src/smt/network_flow_def.h +++ b/src/smt/network_flow_def.h @@ -33,9 +33,7 @@ namespace smt { m_balances.resize(num_nodes); - m_potentials.resize(num_nodes); - m_flows.resize(num_edges); - m_states.resize(num_edges); + m_potentials.resize(num_nodes); m_upwards.resize(num_nodes); m_pred.resize(num_nodes); @@ -58,12 +56,15 @@ namespace smt { m_final[root] = root - 1; m_potentials[root] = numeral::zero(); + m_graph.init_var(root); + fin_numeral sum_supply = fin_numeral::zero(); for (unsigned i = 0; i < m_balances.size(); ++i) { sum_supply += m_balances[i]; } m_balances[root] = -sum_supply; + m_flows.resize(num_nodes + num_edges); m_states.resize(num_nodes + num_edges); m_states.fill(NON_BASIS); @@ -74,14 +75,25 @@ namespace smt { m_depth[i] = 1; m_thread[i] = i + 1; m_final[i] = i; - m_rev_thread[i] = (i == 0) ? root : i - 1; + m_rev_thread[i + 1] = i; 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]; + m_graph.enable_edge(m_graph.add_edge(src, tgt, numeral::one(), explanation())); } TRACE("network_flow", tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread) << pp_vector("Reverse Threads", m_rev_thread) << pp_vector("Last Successors", m_final) << pp_vector("Depths", m_depth);); } + template + bool network_flow::get_edge_id(dl_var source, dl_var target, edge_id & id) { + // m_upwards decides which node is the real source + return m_upwards[source] ? m_graph.get_edge_id(source, target, id) : m_graph.get_edge_id(target, source, id); + } + template void network_flow::update_potentials() { TRACE("network_flow", tout << "update_potentials...\n";); @@ -105,13 +117,13 @@ namespace smt { 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; - m_graph.get_edge_id(u, m_pred[u], e_id); + get_edge_id(u, m_pred[u], e_id); m_flows[e_id] += m_upwards[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; - m_graph.get_edge_id(u, m_pred[u], e_id); + get_edge_id(u, m_pred[u], e_id); m_flows[e_id] += m_upwards[u] ? val : -val; } TRACE("network_flow", tout << pp_vector("Flows", m_flows, true);); @@ -159,6 +171,7 @@ namespace smt { } // Found first common ancestor of source and target m_join_node = u; + TRACE("network_flow", tout << "Found join node " << m_join_node << std::endl;); // FIXME: need to get truly infinite value numeral infty = numeral(INT_MAX); m_delta = infty; @@ -166,7 +179,7 @@ namespace smt { // 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; - m_graph.get_edge_id(u, m_pred[u], e_id); + get_edge_id(u, m_pred[u], e_id); numeral d = m_upwards[u] ? m_flows[e_id] : infty; if (d < m_delta) { m_delta = d; @@ -178,7 +191,7 @@ namespace smt { // 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; - m_graph.get_edge_id(u, m_pred[u], e_id); + get_edge_id(u, m_pred[u], e_id); numeral d = m_upwards[u] ? infty : m_flows[e_id]; if (d <= m_delta) { m_delta = d; @@ -188,8 +201,8 @@ namespace smt { } if (m_delta < infty) { - m_graph.get_edge_id(src, tgt, m_leaving_edge); - TRACE("network_flow", tout << "Found leaving edge" << m_leaving_edge << "between node " << src << " and node " << tgt << "...\n";); + get_edge_id(src, tgt, m_leaving_edge); + TRACE("network_flow", tout << "Found leaving edge " << m_leaving_edge << " between node " << src << " and node " << tgt << "...\n";); return true; } TRACE("network_flow", tout << "Can't find a leaving edge... The problem is unbounded.\n";); @@ -232,8 +245,8 @@ namespace smt { } node gamma = m_thread[m_final[src_in]]; - last = m_pred[gamma] != 0 ? gamma : root; - for (node u = src_in; u == last; u = m_pred[u]) { + last = m_pred[gamma] != -1 ? gamma : root; + for (node u = src_in; u != last; u = m_pred[u]) { m_final[u] = z; } @@ -245,8 +258,8 @@ namespace smt { gamma = m_thread[m_final[tgt_out]]; // REVIEW: check f(u) is not in T_v node delta = m_final[src_out] != m_final[tgt_out] ? m_final[src_out] : m_rev_thread[tgt_out]; - last = m_pred[gamma] != 0 ? gamma : root; - for (node u = src_in; u == last; u = m_pred[u]) { + last = m_pred[gamma] != -1 ? gamma : root; + for (node u = src_in; u != last; u = m_pred[u]) { m_final[u] = delta; } @@ -281,7 +294,7 @@ namespace smt { // Decrease depth of all children in the subtree ++count; int d = m_depth[u] - count; - for (node v = m_thread[u]; v == m_final[u]; v = m_thread[v]) { + for (node v = m_thread[u]; v != m_final[u]; v = m_thread[v]) { m_depth[v] -= d; } } diff --git a/src/smt/theory_diff_logic_def.h b/src/smt/theory_diff_logic_def.h index 70fc681a7..e26947ed5 100644 --- a/src/smt/theory_diff_logic_def.h +++ b/src/smt/theory_diff_logic_def.h @@ -1015,7 +1015,7 @@ bool theory_diff_logic::maximize(theory_var v) { for (unsigned i = 0; i < objective.size(); ++i) { fin_numeral balance(objective[i].second); balances[objective[i].first] = balance; - } + } network_flow net_flow(m_graph, balances); bool is_optimal = net_flow.min_cost();