Add pretty printing for network_flow

Reuse the original graph as much as possible

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(); }

src/smt/network_flow.h

@@ -47,7 +47,7 @@ namespace smt {
         typedef dl_graph<Ext> 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<fin_numeral> m_balances;
@@ -58,9 +58,6 @@ namespace smt {
         // Keep optimal solution of the min cost flow problem
         numeral m_objective_value;
         
-        // Costs on edges
-        vector<fin_numeral> m_costs;
-
         // Basic feasible flows
         vector<numeral> m_flows;
         
@@ -79,18 +76,12 @@ namespace smt {
         svector<node> m_rev_thread;
         // Store a final node of the sub tree rooted at node i
         svector<node> m_final;
-        // Number of nodes in the sub tree rooted at node i
-        svector<int> 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<fin_numeral> const & balances);
-
         // Initialize the network with a feasible spanning tree
         void initialize();
 
@@ -108,12 +99,20 @@ namespace smt {
 
         void update_spanning_tree();
 
-        // Compute the optimal solution
+        std::string pp_vector(std::string const & label, svector<int> v, bool has_header = false);
-        numeral get_optimal_solution(vector<numeral> & result, bool is_dual);
+        std::string pp_vector(std::string const & label, vector<numeral> v, bool has_header = false);
+
+    public:
+
+        network_flow(graph & g, vector<fin_numeral> 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<numeral> & result, bool is_dual);
+
     };
 }
 

src/smt/network_flow_def.h

@@ -28,22 +28,12 @@ namespace smt {
     network_flow<Ext>::network_flow(graph & g, vector<fin_numeral> 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<edge> 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_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<typename Ext>
     void network_flow<Ext>::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_rev_thread[i] = (i == 0) ? root : i - 1;
-            m_num_node[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<typename Ext>
     void network_flow<Ext>::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<typename Ext>
     void network_flow<Ext>::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<typename Ext>
     bool network_flow<Ext>::choose_entering_edge() {
+        TRACE("network_flow", tout << "choose_entering_edge...\n";);
         vector<edge> 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<typename Ext>
     bool network_flow<Ext>::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,7 +159,7 @@ 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; 
         node src, tgt;
@@ -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<typename Ext>
     void network_flow<Ext>::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;
         }
 
-    // Get the optimal solution
+        // Graft T_q and T_r'
-    template<typename Ext>
+        m_thread[x] = src_out;
-    typename network_flow<Ext>::numeral network_flow<Ext>::get_optimal_solution(vector<numeral> & result, bool is_dual) {
+        m_thread[z] = y;
-        m_objective_value = numeral::zero();
+        u = src_out;
-        for (unsigned i = 0; i < m_flows.size(); ++i) {
+        while (u != m_final[src_out]) {
-            m_objective_value += m_costs[i] * m_flows[i];
+            m_depth[u] += 1 + m_depth[src_in];
+            u = m_pred[u];
         }
-        result.reset();
+
-        if (is_dual) {            
+        node gamma = m_thread[m_final[src_in]];
-            result.append(m_potentials);                   
+        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 {
-            result.append(m_flows);     
+                    alpha1 = t1;
+                    alpha2 = t2;
                 }                
-        return m_objective_value;
+                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<typename Ext>
+    std::string network_flow<Ext>::pp_vector(std::string const & label, svector<int> 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>
+    std::string network_flow<Ext>::pp_vector(std::string const & label, vector<numeral> 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();
     }
 
     // Minimize cost flows
@@ -232,6 +345,24 @@ namespace smt {
         }
         return true;
     }
+
+    // Get the optimal solution
+    template<typename Ext>
+    typename network_flow<Ext>::numeral network_flow<Ext>::get_optimal_solution(vector<numeral> & 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

src/smt/theory_diff_logic_def.h

@@ -1008,35 +1008,26 @@ bool theory_diff_logic<Ext>::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<GExt> g;
-    vector<dl_edge<GExt> > const& es = m_graph.get_all_edges();
-    for (unsigned i = 0; i < es.size(); ++i) {
-        dl_edge<GExt> 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<fin_numeral> balances;
+    vector<fin_numeral> 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<GExt> net_flow(g, balances);
+    network_flow<GExt> net_flow(m_graph, balances);
     bool is_optimal = net_flow.min_cost();
     if (is_optimal) {
         vector<numeral> 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;
 }