Add visualization of spanning trees

src/smt/network_flow.h

@@ -85,6 +85,8 @@ namespace smt {
         node m_join_node;
         numeral m_delta;
 
+        unsigned m_step;
+
         // Initialize the network with a feasible spanning tree
         void initialize();
 
@@ -105,6 +107,8 @@ namespace smt {
 
         void update_spanning_tree();
 
+        std::string display_spanning_tree();
+
     public:
 
         network_flow(graph & g, vector<fin_numeral> const & balances);        

src/smt/network_flow_def.h

@@ -59,6 +59,8 @@ namespace smt {
         m_thread.resize(num_nodes);
         m_rev_thread.resize(num_nodes);
         m_final.resize(num_nodes);
+
+        m_step = 0;
     }
 
     template<typename Ext>
@@ -69,6 +71,7 @@ namespace smt {
         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_rev_thread[0] = root;
         m_final[root] = root - 1;
@@ -121,11 +124,12 @@ namespace smt {
                 tout << pp_vector("Depths", m_depth) << pp_vector("Upwards", m_upwards);
                 tout << pp_vector("Potentials", m_potentials) << pp_vector("Flows", m_flows);
             });
+        TRACE("network_flow", tout << "Spanning tree:\n" << display_spanning_tree(););
     }
 
     template<typename Ext>
     edge_id network_flow<Ext>::get_edge_id(dl_var source, dl_var target) {
-        // m_upwards decides which node is the real source
+        // m_upwards[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));
         return id;
@@ -137,8 +141,7 @@ 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]) : 
+        numeral change = m_upwards[src] ? (-cost + m_potentials[src] - m_potentials[tgt]) : (cost - m_potentials[src] + m_potentials[tgt]);
-                            (-cost + m_potentials[src] - m_potentials[tgt]);
         node last = m_thread[m_final[src]];
         for (node u = src; u != last; u = m_thread[u]) {
             m_potentials[u] += change;
@@ -149,7 +152,7 @@ namespace smt {
     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;
+        numeral val = m_states[m_entering_edge] == BASIS ? numeral::zero() : 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]) {
@@ -160,6 +163,7 @@ namespace smt {
         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;
+            TRACE("network_flow", tout << u << ", " << m_pred[u] << ", " << e_id << ", " << m_upwards[u] << ", " << val;);
         }
         TRACE("network_flow", tout << pp_vector("Flows", m_flows, true););
     }
@@ -217,7 +221,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 = get_edge_id(u, m_pred[u]);
-            numeral d = m_upwards[u] ? m_flows[e_id] : infty;
+            numeral d = m_upwards[u] ? infty : m_flows[e_id];
             if (d < m_delta) {
                 m_delta = d;
                 src = u;
@@ -250,75 +254,88 @@ namespace smt {
 
     template<typename Ext>
     void network_flow<Ext>::update_spanning_tree() {        
-        node src_in = m_graph.get_source(m_entering_edge);
+        node p = m_graph.get_source(m_entering_edge);
-        node tgt_in = m_graph.get_target(m_entering_edge);
+        node q = m_graph.get_target(m_entering_edge);
-        node src_out = m_graph.get_source(m_leaving_edge);
+        node u = m_graph.get_source(m_leaving_edge);
-        node tgt_out = m_graph.get_target(m_leaving_edge);
+        node v = m_graph.get_target(m_leaving_edge);
+
         TRACE("network_flow", { 
-            tout << "update_spanning_tree: (" << src_in << ", " << tgt_in << ") enters, (";
+            tout << "update_spanning_tree: (" << p << ", " << q << ") enters, (";
-            tout << src_out << ", " << tgt_out << ") leaves\n";
+            tout << u << ", " << v << ") leaves\n";
         });
-        node root = m_graph.get_num_nodes();
 
-        node x = m_final[src_in];
+        node x = m_final[p];
         node y = m_thread[x];
-        node z = m_final[src_out];
+        node z = m_final[q];
 
-        // Update m_pred (for nodes in the stem from tgt_in to tgt_out)
+        // Update m_pred (for nodes in the stem from q to v)
-        node u = tgt_in;
+        node n = q;
-        node last = m_pred[tgt_out];
+        node last = m_pred[v];
-        node parent = src_in;
+        node parent = p;
-        while (u != last) {
+        while (n != last) {
-            node next = m_pred[u];  
+            node next = m_pred[n];  
-            m_pred[u] = parent;
+            m_pred[n] = parent;
-            m_upwards[u] = !m_upwards[u];
+            m_upwards[n] = !m_upwards[n];
-            parent = u;
+            parent = n;
-            u = next;
+            n = next;
         }
 
+        TRACE("network_flow", tout << "Graft T_q and T_r'\n";);
+
         // Graft T_q and T_r'
-        m_thread[x] = src_out;
+        m_thread[x] = q;
         m_thread[z] = y;
-        u = src_out;
+        n = u;
-        while (u != m_final[src_out]) {
+        last = m_final[u];
-            m_depth[u] += 1 + m_depth[src_in];
+        while (n != last) {
-            u = m_pred[u];
+            m_depth[n] += 1 + m_depth[p];
+            n = m_pred[n];
         }
 
-        node gamma = m_thread[m_final[src_in]];
+        node gamma = m_thread[m_final[p]];        
-        for (node u = src_in; u != gamma; u = m_pred[u]) {
+        n = p;
-            m_final[u] = z;
+        last = m_pred[gamma];
+        while (n != last) {
+            m_final[n] = z;
+            n = m_pred[n];
         }
 
+        TRACE("network_flow", tout << "Update T_r'\n";);
+
         // Update T_r'
-        node phi = m_thread[tgt_out];
+        node phi = m_rev_thread[v];
-        node theta = m_thread[m_final[tgt_out]];
+        node theta = m_thread[m_final[v]];
         m_thread[phi] = theta;
 
-        gamma = m_thread[m_final[tgt_out]];
+        gamma = m_thread[m_final[v]];
-        // REVIEW: check f(u) is not in T_v
+        // REVIEW: check f(n) is not in T_v
-        node delta = m_final[src_out] != m_final[tgt_out] ? m_final[src_out] : m_rev_thread[tgt_out];
+        node delta = m_final[u] != m_final[v] ? m_final[u] : phi;
-        for (node u = src_in; u != gamma; u = m_pred[u]) {
+        n = u;
-            m_final[u] = delta;
+        last = m_pred[gamma];
+        while (n != last) {
+            m_final[n] = delta;
+            n = m_pred[n];
         }
 
         // Reroot T_v at q
-        if (src_out != tgt_in) {
+        if (u != q) {
-            node u = m_pred[src_out];
+            TRACE("network_flow", tout << "Reroot T_v at q\n";);
-            m_thread[m_final[src_out]] = u;
+
-            last = tgt_in;            
+            node n = m_pred[q];
+            m_thread[m_final[q]] = n;
+            last = v;            
             node alpha1, alpha2;
             unsigned count = 0;
-            while (u != last) {
+            while (n != last) {
-                // Find all immediate successors of u
+                // Find all immediate successors of n
-                node t1 = m_thread[u];
+                node t1 = m_thread[n];
                 node t2 = m_thread[m_final[t1]];
                 node t3 = m_thread[m_final[t2]];
-                if (t1 = m_pred[u]) {
+                if (t1 = m_pred[n]) {
                     alpha1 = t2;
                     alpha2 = t3;
                 }
-                else if (t2 = m_pred[u]) {
+                else if (t2 == m_pred[n]) {
                     alpha1 = t1;
                     alpha2 = t3;
                 }
@@ -326,18 +343,18 @@ namespace smt {
                     alpha1 = t1;
                     alpha2 = t2;
                 }                
-                m_thread[u] = alpha1;
+                m_thread[n] = alpha1;
                 m_thread[m_final[alpha1]] = alpha2;
-                u = m_pred[u];
+                n = m_pred[n];
-                m_thread[m_final[alpha2]] = u;
+                m_thread[m_final[alpha2]] = n;
                 // Decrease depth of all children in the subtree
                 ++count;
-                int d = m_depth[u] - count;
+                int d = m_depth[n] - count;
-                for (node v = m_thread[u]; v != m_final[u]; v = m_thread[v]) {
+                for (node m = m_thread[n]; m != m_final[n]; m = m_thread[m]) {
-                    m_depth[v] -= d;
+                    m_depth[m] -= d;
                 }
             }
-            m_thread[m_final[alpha2]] = src_out;
+            m_thread[m_final[alpha2]] = v;
         }
 
         TRACE("network_flow", {
@@ -347,6 +364,36 @@ namespace smt {
             });
     }
 
+    template<typename Ext>
+    std::string network_flow<Ext>::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 << " [" << m_potentials[i] << "]\"];\n";
+        }
+        oss << prefix << root << "[shape=doublecircle,label=\"" << prefix << root << " [" << m_potentials[root] << "]\"];\n";
+        
+        vector<edge> const & es = m_graph.get_all_edges();
+        for (unsigned i = 0; i < es.size(); ++i) {
+            edge const & e = es[i];
+            if (e.is_enabled()) {
+                oss << prefix << e.get_source() << " -> " << prefix << e.get_target();
+                if (m_states[i] == BASIS) {
+                    oss << "[color=red,penwidth=3.0,label=\"" << m_flows[i] << "\"];\n";
+                }
+                else {
+                    oss << "[label=\"" << m_flows[i] << "\"];\n";
+                }
+            }
+        }
+        oss << std::endl;
+        return oss.str();
+    }
+
     // Minimize cost flows
     // Return true if found an optimal solution, and return false if unbounded
     template<typename Ext>
@@ -355,11 +402,13 @@ namespace smt {
         while (choose_entering_edge()) {
             bool bounded = choose_leaving_edge();
             if (!bounded) return false;
+            update_flows();
             if (m_entering_edge != m_leaving_edge) {
                 m_states[m_entering_edge] = BASIS;
                 m_states[m_leaving_edge] = NON_BASIS;
                 update_spanning_tree();
                 update_potentials();                
+                TRACE("network_flow", tout << "Spanning tree:\n" << display_spanning_tree(););
             }
         }
         TRACE("network_flow", tout << "Found optimal solution.\n";);

src/smt/theory_diff_logic_def.h

@@ -1007,7 +1007,7 @@ bool theory_diff_logic<Ext>::maximize(theory_var v) {
         for (unsigned i = 0; i < objective.size(); ++i) {
             verbose_stream() << "coefficient " << objective[i].second << " of theory_var " << objective[i].first << "\n";
         }
-        verbose_stream() << m_objective_consts[v] << "\n";);
+        verbose_stream() << "free coefficient " << m_objective_consts[v] << "\n";);
 
     // Objective coefficients now become balances
     vector<fin_numeral> balances(m_graph.get_num_nodes());
@@ -1024,9 +1024,9 @@ bool theory_diff_logic<Ext>::maximize(theory_var v) {
         m_objective_value = net_flow.get_optimal_solution(potentials, true);
         std::cout << "Objective value of var " << v << ": " << m_objective_value << std::endl;
         // TODO: return the model of the optimal solution from potentials
-        IF_VERBOSE(1,
+        TRACE("network_flow",
             for (unsigned i = 0; i < potentials.size(); ++i) {
-                verbose_stream() << "v" << i << " -> " << potentials[i] << "\n";
+                tout << "v" << i << " -> " << potentials[i] << "\n";
             });
 
     }