Reduce difference logic solver to min cost flow

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.
+    // If there is such edge, return its edge_id in parameter id.
-    bool get_edge(dl_var source, dl_var target, edge & e) {
+    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];
+            edge & e    = m_edges[e_id];
-            if (e0.is_enabled() && e0.get_target() == target && !found) {
+            if (e.is_enabled() && e.get_target() == target && !found) {
-                e     = e0;
+                id    = e_id;
                 found = true;
             }
         }
         return found;
     }
 
+    edges & get_all_edges() {
+        return m_edges;
+    }
+
     template<typename Functor>
     void enumerate_edges(dl_var source, dl_var target, Functor& f) {
         edge_id_vector & edges = m_out_edges[source];

src/smt/network_flow.h

@@ -38,31 +38,29 @@ namespace smt {
     // Solve minimum cost flow problem using Network Simplex algorithm
     template<typename Ext>
     class network_flow : private Ext {
-        struct GExt : public Ext {
-            typedef literal explanation;
-        };
-
         typedef dl_var node;
-        typedef dl_edge<GExt> edge;
+        typedef dl_edge<Ext> edge;
-        typedef dl_graph<GExt> graph;     
+        typedef dl_graph<Ext> graph;     
-        typedef typename Ext::numeral numeral;        
+        typedef typename Ext::numeral numeral;
         graph m_graph;
 
         // Denote supply/demand b_i on node i
         vector<numeral> m_balances;
 
+        // Duals of flows which are convenient to compute dual solutions
         vector<numeral> m_potentials;
 
         // Keep optimal solution of the min cost flow problem
         inf_rational m_objective;
 
+        // Costs on edges
+        vector<numeral> & m_costs;
+
         // Basic feasible flows
         vector<numeral> m_flows;
 
-        // Denote the spanning tree of basic edges
+        // An element is true if the corresponding edge points upwards (compared to the root node)
-        vector<edge> m_basics;
+        svector<bool> m_upwards;
-        // Denote non-basic edges with flow 0 for uncapicitated networks
-        vector<edge> m_nonbasics;
 
         // Store the parent of a node in the spanning tree
         svector<node> m_pred;
@@ -71,7 +69,12 @@ namespace smt {
         // Store the pointer to the next node in depth first search ordering
         svector<node> m_thread;
 
+        bool m_is_optimal;
+
     public:
+
+        network_flow(graph & g, vector<numeral> & 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);
+        edge_id choose_leaving_edge(edge_id entering_edge);
-
-        void update_basics(const edge & entering_edge, const edge & leaving_edge);
 
+        void update_spanning_tree(edge_id entering_edge, edge_id leaving_edge);
+        
         bool is_unbounded();
 
         // Compute the optimal solution

src/smt/network_flow_def.h

@@ -24,6 +24,12 @@ Notes:
 
 namespace smt {
 
+    template<typename Ext>
+    network_flow<Ext>::network_flow(graph & g, vector<numeral> & costs) :
+        m_graph(g),
+        m_costs(costs) {
+    }
+
     template<typename Ext>
     void network_flow<Ext>::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<rational, m_potentials.size()> potentials;
-        std::copy(m_potentials.begin(), m_potentials.end(), potentials);
-        rational zero(0);
-        potentials[0] = zero;
-        node next = m_thread[0];
         
-        while (next != 0) {
+        numeral zero(0);
-            node current = m_pred[next];
+        m_potentials.set(0, zero);
-            edge e;
+        node target = m_thread[0];
-            if (m_graph.get_edge(current, next, e)) {
+        
-                potentials[next] = potentials[current] - e.get_weight();
+        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)) {
+            if (m_graph.get_edge_id(target, source, e_id)) {
-                potentials[next] = potentials[current] + e.get_weight();
+                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<typename Ext>
     void network_flow<Ext>::compute_flows() {
         vector<numeral> balances(m_balances);
-        numeral zero;
+        
-        m_flows.fill(zero);
+        // OPTIMIZE: Need a set data structure for efficiently removing elements
-        vector<edge> basics(m_basics);
+        vector<edge_id> basics;
-        // TODO: need a way to find a leaf node of a spanning tree
         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<typename Ext>
-    bool network_flow<Ext>::is_optimal(edge & violating_edge) {        
+    bool network_flow<Ext>::is_optimal(edge_id & violating_edge) {    
-        typename vector<edge>::iterator it  = m_nonbasics.begin();
+        // TODO: how to get nonbasics vector?
-        typename vector<edge>::iterator end = m_nonbasics.end();
+        vector<edge> nonbasics;
+        typename vector<edge>::iterator it  = nonbasics.begin();
+        typename vector<edge>::iterator end = nonbasics.end();
         bool found = false;
-        for (unsigned int i = 0; i < m_nonbasics.size(); ++i) {
+        for (unsigned int i = 0; i < nonbasics.size(); ++i) {
-            edge & e = m_nonbasics[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) {
+                numeral zero(0);
-                    violating_edge = e;
+                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<typename Ext>
-    dl_edge<typename network_flow<Ext>::GExt> network_flow<Ext>::choose_leaving_edge(const edge & entering_edge) {
+    edge_id network_flow<Ext>::choose_leaving_edge(edge_id entering_edge) {
-        node source = entering_edge.get_source();
+        node source = m_graph.get_source(entering_edge);
-        node target = entering_edge.get_target();
+        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;
+        edge_id e_id;
-        m_graph.get_edge(source, target, e);
+        m_graph.get_edge_id(source, target, e_id);
-        return e;
+        return e_id;
     }
 
     template<typename Ext>
-    void network_flow<Ext>::update_basics(const edge & entering_edge, const edge & leaving_edge) {
+    void network_flow<Ext>::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<typename Ext>
@@ -124,24 +164,34 @@ namespace smt {
     // Get the optimal solution
     template<typename Ext>
     void network_flow<Ext>::get_optimal_solution(numeral & objective, vector<numeral> & 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<typename Ext>
     bool network_flow<Ext>::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();
+            edge_id leaving_edge = choose_leaving_edge(entering_edge);
-            update_tree(entering_edge, leaving_edge);
+            update_spanning_tree(entering_edge, leaving_edge);
-            if (is_unbounded())
+            if (is_unbounded()) {
-                return false;
+                m_is_optimal = false;
+                return m_is_optimal;
+            }
         }
-        return true;
+        m_is_optimal = true;
+        return m_is_optimal;
     }
 }
 

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<inf_rational> get_objective_value(theory_var v);
+        numeral m_objective_value;
 
         typedef vector <std::pair<theory_var, rational> > objective_term;
         vector<objective_term> m_objectives;
 
         void internalize_objective(app * n, objective_term & objective);
 
-        network_flow<Ext> m_network_flow;
-
     private:        
 
         virtual void new_eq_eh(theory_var v1, theory_var v2, justification& j);

src/smt/theory_diff_logic_def.h

@@ -1005,7 +1005,37 @@ bool theory_diff_logic<Ext>::maximize(theory_var v) {
                 }
                 verbose_stream() << "\n";);
     NOT_IMPLEMENTED_YET();
-    return false;
+    // Double the number of edges in the new graph
+    dl_graph<GExt> g;
+    vector<dl_edge<GExt>> es = m_graph.get_all_edges();
+    dl_var offset = m_graph.get_num_edges();
+    for (unsigned i = 0; i < es.size(); ++i) {
+        dl_edge<GExt> 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<numeral> 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<numeral> costs;
+    costs.append(base_costs);
+    costs.append(aux_costs);
+    
+    network_flow<GExt> net_flow(g, costs);
+    bool is_optimal = net_flow.min_cost();
+    if (is_optimal) {
+        numeral objective_value;
+        vector<numeral> 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<typename Ext>
@@ -1018,6 +1048,7 @@ theory_var theory_diff_logic<Ext>::add_objective(app* term) {
 
 template<typename Ext>
 inf_eps_rational<inf_rational> theory_diff_logic<Ext>::get_objective_value(theory_var v) { 
+    NOT_IMPLEMENTED_YET();
 	inf_rational objective;
     inf_eps_rational<inf_rational> val(objective);     
     return val;