From 532c345fd119e0b2b08534c171cd0ca568040056 Mon Sep 17 00:00:00 2001 From: Anh-Dung Phan Date: Fri, 25 Oct 2013 17:42:03 -0700 Subject: [PATCH 1/2] Reduce difference logic solver to min cost flow --- src/smt/diff_logic.h | 14 ++-- src/smt/network_flow.h | 33 ++++---- src/smt/network_flow_def.h | 134 ++++++++++++++++++++++---------- src/smt/theory_diff_logic.h | 3 +- src/smt/theory_diff_logic_def.h | 33 +++++++- 5 files changed, 152 insertions(+), 65 deletions(-) diff --git a/src/smt/diff_logic.h b/src/smt/diff_logic.h index 74a5dbbd2..23f0de06c 100644 --- a/src/smt/diff_logic.h +++ b/src/smt/diff_logic.h @@ -932,23 +932,27 @@ public: } // Return true if there is an edge source --> target. - // If there is such edge, return it in parameter e. - bool get_edge(dl_var source, dl_var target, edge & e) { + // If there is such edge, return its edge_id in parameter id. + bool get_edge_id(dl_var source, dl_var target, edge_id & id) { edge_id_vector & edges = m_out_edges[source]; typename edge_id_vector::iterator it = edges.begin(); typename edge_id_vector::iterator end = edges.end(); bool found = false; for (; it != end; ++it) { edge_id e_id = *it; - edge & e0 = m_edges[e_id]; - if (e0.is_enabled() && e0.get_target() == target && !found) { - e = e0; + edge & e = m_edges[e_id]; + if (e.is_enabled() && e.get_target() == target && !found) { + id = e_id; found = true; } } return found; } + edges & get_all_edges() { + return m_edges; + } + template void enumerate_edges(dl_var source, dl_var target, Functor& f) { edge_id_vector & edges = m_out_edges[source]; diff --git a/src/smt/network_flow.h b/src/smt/network_flow.h index d10ccf005..8232f4d0a 100644 --- a/src/smt/network_flow.h +++ b/src/smt/network_flow.h @@ -38,31 +38,29 @@ namespace smt { // Solve minimum cost flow problem using Network Simplex algorithm template class network_flow : private Ext { - struct GExt : public Ext { - typedef literal explanation; - }; - typedef dl_var node; - typedef dl_edge edge; - typedef dl_graph graph; - typedef typename Ext::numeral numeral; + typedef dl_edge edge; + typedef dl_graph graph; + typedef typename Ext::numeral numeral; graph m_graph; // Denote supply/demand b_i on node i vector m_balances; + // Duals of flows which are convenient to compute dual solutions vector m_potentials; // Keep optimal solution of the min cost flow problem inf_rational m_objective; + // Costs on edges + vector & m_costs; + // Basic feasible flows vector m_flows; - // Denote the spanning tree of basic edges - vector m_basics; - // Denote non-basic edges with flow 0 for uncapicitated networks - vector m_nonbasics; + // An element is true if the corresponding edge points upwards (compared to the root node) + svector m_upwards; // Store the parent of a node in the spanning tree svector m_pred; @@ -71,7 +69,12 @@ namespace smt { // Store the pointer to the next node in depth first search ordering svector m_thread; + bool m_is_optimal; + public: + + network_flow(graph & g, vector & costs); + // Initialize the network with a feasible spanning tree void initialize(); @@ -81,13 +84,13 @@ namespace smt { // If all reduced costs are non-negative, the current flow is optimal // If not optimal, return a violating edge in the corresponding variable - bool is_optimal(edge & violating_edge); + bool is_optimal(edge_id & violating_edge); // Send as much flow as possible around the cycle, the first basic edge with flow 0 will leave - edge choose_leaving_edge(const edge & entering_edge); - - void update_basics(const edge & entering_edge, const edge & leaving_edge); + edge_id choose_leaving_edge(edge_id entering_edge); + void update_spanning_tree(edge_id entering_edge, edge_id leaving_edge); + bool is_unbounded(); // Compute the optimal solution diff --git a/src/smt/network_flow_def.h b/src/smt/network_flow_def.h index 8ce0a2458..3c5106ac6 100644 --- a/src/smt/network_flow_def.h +++ b/src/smt/network_flow_def.h @@ -24,6 +24,12 @@ Notes: namespace smt { + template + network_flow::network_flow(graph & g, vector & costs) : + m_graph(g), + m_costs(costs) { + } + template void network_flow::initialize() { // TODO: construct an initial spanning tree i.e. inializing m_pred, m_depth and m_thread. @@ -36,54 +42,74 @@ namespace smt { SASSERT(!m_potentials.empty()); SASSERT(!m_thread.empty()); SASSERT(m_thread.size() == m_pred.size()); - - array potentials; - std::copy(m_potentials.begin(), m_potentials.end(), potentials); - rational zero(0); - potentials[0] = zero; - node next = m_thread[0]; - while (next != 0) { - node current = m_pred[next]; - edge e; - if (m_graph.get_edge(current, next, e)) { - potentials[next] = potentials[current] - e.get_weight(); + numeral zero(0); + m_potentials.set(0, zero); + node target = m_thread[0]; + + while (target != 0) { + node source = m_pred[target]; + edge_id e_id; + if (m_graph.get_edge_id(source, target, e_id)) { + m_potentials.set(target, m_potentials[source] - m_graph.get_weight(e_id)); } - if (m_graph.get_edge(next, current, e)) { - potentials[next] = potentials[current] + e.get_weight(); + if (m_graph.get_edge_id(target, source, e_id)) { + m_potentials.set(target, m_potentials[source] + m_graph.get_weight(e_id)); } - next = m_thread[next]; + target = m_thread[target]; } - std::copy(potentials.begin(), potentials.end(), m_potentials); } template void network_flow::compute_flows() { vector balances(m_balances); - numeral zero; - m_flows.fill(zero); - vector basics(m_basics); - // TODO: need a way to find a leaf node of a spanning tree + + // OPTIMIZE: Need a set data structure for efficiently removing elements + vector basics; while (!basics.empty()) { - return; + // Find a leaf node of a spanning tree + node target; + for (unsigned int i = 0; i < m_thread.size(); ++i) { + if (m_depth[i] <= m_depth[m_thread[i]]) { + target = i; + break; + } + } + node source = m_pred[target]; + edge_id e_id; + if (m_graph.get_edge_id(source, target, e_id)) { + m_flows.set(e_id, -m_graph.get_weight(basics[target])); + basics[source] += basics[target]; + basics.erase(e_id); + } + else if (m_graph.get_edge_id(target, source, e_id)) { + m_flows.set(e_id, m_graph.get_weight(basics[target])); + basics[source] += basics[target]; + basics.erase(e_id); + } } } template - bool network_flow::is_optimal(edge & violating_edge) { - typename vector::iterator it = m_nonbasics.begin(); - typename vector::iterator end = m_nonbasics.end(); + bool network_flow::is_optimal(edge_id & violating_edge) { + // TODO: how to get nonbasics vector? + vector nonbasics; + typename vector::iterator it = nonbasics.begin(); + typename vector::iterator end = nonbasics.end(); bool found = false; - for (unsigned int i = 0; i < m_nonbasics.size(); ++i) { - edge & e = m_nonbasics[i]; + for (unsigned int i = 0; i < nonbasics.size(); ++i) { + edge & e = nonbasics[i]; if (e.is_enabled()) { node source = e.get_source(); node target = e.get_target(); numeral cost = e.get_weight() - m_potentials[source] + m_potentials[target]; // Choose the first negative-cost edge to be the violating edge // TODO: add multiple pivoting strategies - if (cost < 0) { - violating_edge = e; + numeral zero(0); + if (cost < zero) { + edge_id e_id; + m_graph.get_edge_id(source, target, e_id); + violating_edge = e_id; found = true; break; } @@ -93,9 +119,9 @@ namespace smt { } template - dl_edge::GExt> network_flow::choose_leaving_edge(const edge & entering_edge) { - node source = entering_edge.get_source(); - node target = entering_edge.get_target(); + edge_id network_flow::choose_leaving_edge(edge_id entering_edge) { + node source = m_graph.get_source(entering_edge); + node target = m_graph.get_target(entering_edge); while (source != target) { if (m_depth[source] > m_depth[target]) source = m_pred[source]; @@ -106,14 +132,28 @@ namespace smt { target = m_pred[target]; } } - edge e; - m_graph.get_edge(source, target, e); - return e; + edge_id e_id; + m_graph.get_edge_id(source, target, e_id); + return e_id; } template - void network_flow::update_basics(const edge & entering_edge, const edge & leaving_edge) { - + void network_flow::update_spanning_tree(edge_id entering_edge, edge_id leaving_edge) { + // Need special handling in case two edges are identical + SASSERT(entering_edge != leaving_edge); + + // Update potentials + node target = m_upwards[leaving_edge] ? m_graph.get_source(leaving_edge) : m_graph.get_target(leaving_edge); + numeral src_pot = m_potentials[m_graph.get_source(entering_edge)]; + numeral tgt_pot = m_potentials[m_graph.get_target(entering_edge)]; + numeral weight = m_graph.get_weight(entering_edge); + numeral change = m_upwards[entering_edge] ? (weight - src_pot + tgt_pot) : (-weight + src_pot - tgt_pot); + m_potentials[target] += change; + node start = m_thread[target]; + while (m_depth[start] > m_depth[target]) { + m_potentials[start] += change; + start = m_thread[start]; + } } template @@ -124,24 +164,34 @@ namespace smt { // Get the optimal solution template void network_flow::get_optimal_solution(numeral & objective, vector & flows) { + SASSERT(m_is_optimal); flows.reset(); flows.append(m_flows); - // TODO: calculate objective value + numeral cost(0); + for (unsigned int i = 0; i < m_flows.size(); ++i) { + // FIXME: this * operator is not supported + //cost += m_costs[i] * m_flows[i]; + } + objective = cost; } // Minimize cost flows // Return true if found an optimal solution, and return false if unbounded template bool network_flow::min_cost() { + SASSERT(!m_graph.get_all_edges().empty()); initialize(); - edge & entering_edge; + edge_id entering_edge; while (!is_optimal(entering_edge)) { - edge & leaving_edge = choose_leaving_edge(); - update_tree(entering_edge, leaving_edge); - if (is_unbounded()) - return false; + edge_id leaving_edge = choose_leaving_edge(entering_edge); + update_spanning_tree(entering_edge, leaving_edge); + if (is_unbounded()) { + m_is_optimal = false; + return m_is_optimal; + } } - return true; + m_is_optimal = true; + return m_is_optimal; } } diff --git a/src/smt/theory_diff_logic.h b/src/smt/theory_diff_logic.h index 9447b0d3e..3af150e64 100644 --- a/src/smt/theory_diff_logic.h +++ b/src/smt/theory_diff_logic.h @@ -307,14 +307,13 @@ namespace smt { virtual bool maximize(theory_var v); virtual theory_var add_objective(app* term); virtual inf_eps_rational get_objective_value(theory_var v); + numeral m_objective_value; typedef vector > objective_term; vector m_objectives; void internalize_objective(app * n, objective_term & objective); - network_flow m_network_flow; - private: virtual void new_eq_eh(theory_var v1, theory_var v2, justification& j); diff --git a/src/smt/theory_diff_logic_def.h b/src/smt/theory_diff_logic_def.h index 684850d77..841375253 100644 --- a/src/smt/theory_diff_logic_def.h +++ b/src/smt/theory_diff_logic_def.h @@ -1005,7 +1005,37 @@ bool theory_diff_logic::maximize(theory_var v) { } verbose_stream() << "\n";); NOT_IMPLEMENTED_YET(); - return false; + // Double the number of edges in the new graph + dl_graph g; + vector> es = m_graph.get_all_edges(); + dl_var offset = m_graph.get_num_edges(); + for (unsigned i = 0; i < es.size(); ++i) { + dl_edge e(es[i]); + g.enable_edge(g.add_edge(e)); + g.enable_edge(g.add_edge(e.get_target() + offset, e.get_source() + offset, e.get_weight(), e.get_explanation())); + } + + // Objective coefficients now become costs + vector base_costs, aux_costs; + for (unsigned i = 0; i < m_objectives[v].size(); ++i) { + numeral cost(m_objectives[v][i].second); + base_costs.push_back(cost); + aux_costs.push_back(-cost); + } + vector costs; + costs.append(base_costs); + costs.append(aux_costs); + + network_flow net_flow(g, costs); + bool is_optimal = net_flow.min_cost(); + if (is_optimal) { + numeral objective_value; + vector flows; + net_flow.get_optimal_solution(objective_value, flows); + m_objective_value = objective_value.get_rational(); + // TODO: return the model of the optimal solution + } + return is_optimal; } template @@ -1018,6 +1048,7 @@ theory_var theory_diff_logic::add_objective(app* term) { template inf_eps_rational theory_diff_logic::get_objective_value(theory_var v) { + NOT_IMPLEMENTED_YET(); inf_rational objective; inf_eps_rational val(objective); return val; From 3d943bf70dc40ed8e02f185ace60c53c8dcfecc6 Mon Sep 17 00:00:00 2001 From: Anh-Dung Phan Date: Sat, 26 Oct 2013 05:22:52 +0200 Subject: [PATCH 2/2] Fix a mistake in previous commit causing imcompilable code Also correct my alias --- src/smt/theory_diff_logic_def.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/smt/theory_diff_logic_def.h b/src/smt/theory_diff_logic_def.h index 841375253..41fbf7566 100644 --- a/src/smt/theory_diff_logic_def.h +++ b/src/smt/theory_diff_logic_def.h @@ -1010,8 +1010,8 @@ bool theory_diff_logic::maximize(theory_var v) { vector> es = m_graph.get_all_edges(); dl_var offset = m_graph.get_num_edges(); for (unsigned i = 0; i < es.size(); ++i) { - dl_edge e(es[i]); - g.enable_edge(g.add_edge(e)); + dl_edge & e = es[i]; + g.enable_edge(g.add_edge(e.get_source(), e.get_target(), e.get_weight(), e.get_explanation())); g.enable_edge(g.add_edge(e.get_target() + offset, e.get_source() + offset, e.get_weight(), e.get_explanation())); }