Add a vector of edges to handle spanning trees

src/smt/network_flow.h

@@ -51,7 +51,7 @@ namespace smt {
         typedef typename Ext::fin_numeral fin_numeral;
 
         graph m_graph;
-        spanning_tree_base m_tree;
+        thread_spanning_tree<Ext> m_tree;
 
         // Denote supply/demand b_i on node i
         vector<fin_numeral> m_balances;
@@ -68,7 +68,6 @@ namespace smt {
 
         edge_id m_enter_id, m_leave_id;
         optional<numeral> m_delta;
-        bool m_is_swap_enter, m_is_swap_leave;
 
         // Initialize the network with a feasible spanning tree
         void initialize();

src/smt/network_flow_def.h

@@ -27,7 +27,8 @@ namespace smt {
 
     template<typename Ext>
     network_flow<Ext>::network_flow(graph & g, vector<fin_numeral> const & balances) :
-        m_balances(balances) {
+        m_balances(balances),
+        m_tree(m_graph) {
         // Network flow graph has the edges in the reversed order compared to constraint graph
         // We only take enabled edges from the original graph
         for (unsigned i = 0; i < g.get_num_nodes(); ++i) {
@@ -40,14 +41,13 @@ namespace smt {
                 m_graph.add_edge(e.get_target(), e.get_source(), e.get_weight(), explanation());
             }
         }
-        m_tree = thread_spanning_tree<Ext>(m_graph);
         m_step = 0;
     }
 
     template<typename Ext>
     void network_flow<Ext>::initialize() {
         TRACE("network_flow", tout << "initialize...\n";);
-        // Create an artificial root node to construct initial spanning m_tree
+        // Create an artificial root node to construct initial spanning tree
         unsigned num_nodes = m_graph.get_num_nodes();
         unsigned num_edges = m_graph.get_num_edges();
         
@@ -69,24 +69,25 @@ namespace smt {
         m_states.resize(num_nodes + num_edges);
         m_states.fill(LOWER);
 
-        // Create artificial edges from/to root node to/from other nodes and initialize the spanning m_tree
+        // Create artificial edges from/to root node to/from other nodes and initialize the spanning tree
-        svector<bool> upwards(num_nodes, false);
+        svector<edge_id> tree;
+        bool is_forward;
         for (unsigned i = 0; i < num_nodes; ++i) {
-            upwards[i] = !m_balances[i].is_neg();
+            is_forward = !m_balances[i].is_neg();
             m_states[num_edges + i] = BASIS;
-            node src = upwards[i] ? i : root;
+            node src = is_forward ? i : root;
-            node tgt = upwards[i] ? root : i;            
+            node tgt = is_forward ? root : i;            
-            m_flows[num_edges + i] = upwards[i] ? m_balances[i] : -m_balances[i];
+            m_flows[num_edges + i] = is_forward ? m_balances[i] : -m_balances[i];
-            m_potentials[i] = upwards[i] ? numeral::one() : -numeral::one();
+            m_potentials[i] = is_forward ? numeral::one() : -numeral::one();
-            m_graph.add_edge(src, tgt, numeral::one(), explanation());
+            tree.push_back(m_graph.add_edge(src, tgt, numeral::one(), explanation()));
         }
 
-        m_tree.initialize(upwards);
+        m_tree.initialize(tree);
 
         TRACE("network_flow", {
                 tout << pp_vector("Potentials", m_potentials, true) << pp_vector("Flows", m_flows);
             });
-        TRACE("network_flow", tout << "Spanning m_tree:\n" << display_spanning_tree(););
+        TRACE("network_flow", tout << "Spanning tree:\n" << display_spanning_tree(););
         SASSERT(check_well_formed());
     }
 
@@ -95,11 +96,15 @@ namespace smt {
         node src = m_graph.get_source(m_enter_id);
         node tgt = m_graph.get_target(m_enter_id);      
         numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(m_enter_id);
-        numeral change = m_is_swap_leave ? -cost : cost;
+        numeral change = m_tree.in_subtree_t2(tgt) ? cost : -cost;
+        node start = m_graph.get_source(m_leave_id);
+        if (!m_tree.in_subtree_t2(start)) {
+            start = m_graph.get_target(m_leave_id);;
+        }
+        TRACE("network_flow", tout << "update_potentials of T_" << start << " with change = " << change << "...\n";);
         svector<node> descendants;
-        node start = m_is_swap_enter ? src : tgt;
-        TRACE("network_flow", tout << "update_potentials of T_" << start << " with delta = " << change << "...\n";);
         m_tree.get_descendants(start, descendants);
+        SASSERT(descendants.size() >= 1);
         for (unsigned i = 0; i < descendants.size(); ++i) {
             node u = descendants[i];           
             m_potentials[u] += change;
@@ -110,16 +115,16 @@ namespace smt {
     template<typename Ext>
     void network_flow<Ext>::update_flows() {
         TRACE("network_flow", tout << "update_flows...\n";);
-        numeral val = *m_delta;
+        m_flows[m_enter_id] += *m_delta;
-        m_flows[m_enter_id] += val;
         node src = m_graph.get_source(m_enter_id);
         node tgt = m_graph.get_target(m_enter_id); 
         svector<edge_id> path;
         svector<bool> against;
         m_tree.get_path(src, tgt, path, against);
+        SASSERT(path.size() >= 1);
         for (unsigned i = 0; i < path.size(); ++i) {
             edge_id e_id = path[i];
-            m_flows[e_id] += against[i] ? -val : val;
+            m_flows[e_id] += against[i] ? - *m_delta : *m_delta;
         }
         TRACE("network_flow", tout << pp_vector("Flows", m_flows, true););
     }
@@ -127,8 +132,8 @@ namespace smt {
     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();
+        unsigned num_edges = m_graph.get_num_edges();
-        for (unsigned i = 0; i < es.size(); ++i) {
+        for (unsigned i = 0; i < num_edges; ++i) {
             node src = m_graph.get_source(i);
             node tgt = m_graph.get_target(i);
             if (m_states[i] != BASIS) {
@@ -138,7 +143,7 @@ namespace smt {
                     m_enter_id = i;
                     TRACE("network_flow", {
                         tout << "Found entering edge " << i << " between node ";
-                        tout << src << " and node " << tgt << "...\n";
+                        tout << src << " and node " << tgt << " with reduced cost = " << cost << "...\n";
                     });
                     return true;
                 }
@@ -158,9 +163,10 @@ namespace smt {
         svector<edge_id> path;
         svector<bool> against;
         m_tree.get_path(src, tgt, path, against);
+        SASSERT(path.size() >= 1);
         for (unsigned i = 0; i < path.size(); ++i) {
             edge_id e_id = path[i];
-            if (against[i] && (!m_delta || m_flows[e_id] <= *m_delta)) {
+            if (against[i] && (!m_delta || m_flows[e_id] < *m_delta)) {
                 m_delta = m_flows[e_id];
                 leave_id = e_id;
             }
@@ -182,7 +188,7 @@ namespace smt {
 
     template<typename Ext>
     void network_flow<Ext>::update_spanning_tree() {  
-        m_tree.update(m_enter_id, m_leave_id, m_is_swap_enter, m_is_swap_leave);
+        m_tree.update(m_enter_id, m_leave_id);
     }
 
     // Minimize cost flows
@@ -201,7 +207,7 @@ namespace smt {
                 m_states[m_leave_id] = (m_flows[m_leave_id].is_zero()) ? LOWER : UPPER;
                 update_spanning_tree();
                 update_potentials();                
-                TRACE("network_flow", tout << "Spanning m_tree:\n" << display_spanning_tree(););
+                TRACE("network_flow", tout << "Spanning tree:\n" << display_spanning_tree(););
                 SASSERT(check_well_formed());
             } 
             else {
@@ -217,12 +223,11 @@ namespace smt {
     template<typename Ext>
     typename network_flow<Ext>::numeral network_flow<Ext>::get_optimal_solution(vector<numeral> & result, bool is_dual) {
         numeral objective_value = numeral::zero();
-        vector<edge> const & es = m_graph.get_all_edges();
+        unsigned num_edges = m_graph.get_num_edges();
-        for (unsigned i = 0; i < es.size(); ++i) {
+        for (unsigned i = 0; i < num_edges; ++i) {
-            edge const & e = es[i];
             if (m_states[i] == BASIS) 
             {
-                objective_value += e.get_weight().get_rational() * m_flows[i];
+                objective_value += m_graph.get_weight(i).get_rational() * m_flows[i];
             }
         }
         result.reset();
@@ -243,6 +248,7 @@ namespace smt {
     template<typename Ext>
     bool network_flow<Ext>::check_well_formed() {
         SASSERT(m_tree.check_well_formed());
+        SASSERT(!m_delta || !(*m_delta).is_neg());
 
         // m_flows are zero on non-basic edges
         for (unsigned i = 0; i < m_flows.size(); ++i) {
@@ -250,11 +256,10 @@ namespace smt {
             SASSERT(m_states[i] == BASIS || m_flows[i].is_zero());
         }
 
-        vector<edge> const & es = m_graph.get_all_edges();
+        unsigned num_edges = m_graph.get_num_edges();
-        for (unsigned i = 0; i < es.size(); ++i) {
+        for (unsigned i = 0; i < num_edges; ++i) {
-            edge const & e = es[i];
             if (m_states[i] == BASIS) {
-                SASSERT(m_potentials[e.get_source()] - m_potentials[e.get_target()] == e.get_weight());
+                SASSERT(m_potentials[m_graph.get_source(i)] - m_potentials[m_graph.get_target(i)] == m_graph.get_weight(i));
             }
         }
 
@@ -264,11 +269,10 @@ namespace smt {
     template<typename Ext>
     bool network_flow<Ext>::check_optimal() {
         numeral total_cost = numeral::zero();
-        vector<edge> const & es = m_graph.get_all_edges();
+        unsigned num_edges = m_graph.get_num_edges();
-        for (unsigned i = 0; i < es.size(); ++i) {
+        for (unsigned i = 0; i < num_edges; ++i) {
-            edge const & e = es[i];
             if (m_states[i] == BASIS) {
-                total_cost += e.get_weight().get_rational() * m_flows[i];
+                total_cost += m_graph.get_weight(i).get_rational() * m_flows[i];
             }
         }
 
@@ -299,15 +303,14 @@ namespace smt {
             oss << prefix << root << "[shape=doublecircle,label=\"" << prefix << root << " [";
             oss << m_potentials[root] << "/" << m_balances[root] << "]\"];\n";
         
-            vector<edge> const & es = m_graph.get_all_edges();
+            unsigned num_edges = m_graph.get_num_edges();
-            for (unsigned i = 0; i < es.size(); ++i) {
+            for (unsigned i = 0; i < num_edges; ++i) {
-                edge const & e = es[i];
+                oss << prefix << m_graph.get_source(i) << " -> " << prefix << m_graph.get_target(i);
-                oss << prefix << e.get_source() << " -> " << prefix << e.get_target();
                 if (m_states[i] == BASIS) {
-                    oss << "[color=red,penwidth=3.0,label=\"" << m_flows[i] << "/" << e.get_weight() << "\"];\n";
+                    oss << "[color=red,penwidth=3.0,label=\"" << m_flows[i] << "/" << m_graph.get_weight(i) << "\"];\n";
                 }
                 else {
-                    oss << "[label=\"" << m_flows[i] << "/" << e.get_weight() << "\"];\n";
+                    oss << "[label=\"" << m_flows[i] << "/" << m_graph.get_weight(i) << "\"];\n";
                 }
             }
             oss << std::endl;

src/smt/spanning_tree.h

@@ -40,9 +40,10 @@ namespace smt {
         // Store the pointer from node i to the next node in depth-first search order
         svector<node> m_thread;
 
-        // m_upwards[i] is true if the corresponding edge 
+        // i |-> edge between (i, m_pred[i])
-        // (i, m_pred[i]) points upwards (pointing toward the root node)
+        svector<edge_id> m_tree;
-        svector<bool> m_upwards;
+
+        node m_root_t2;
 
         graph & m_graph;
 
@@ -51,20 +52,23 @@ namespace smt {
         void fix_depth(node start, node end);
         node get_final(int start);
         bool is_preorder_traversal(node start, node end);   
-        edge_id get_edge_to_parent(node start) const;
         node get_common_ancestor(node u, node v);
-
-    public:
-        thread_spanning_tree(graph & g);
-
-        void initialize(svector<bool> const & upwards);
-        void get_descendants(node start, svector<node> & descendants);
-        
-        void update(edge_id enter_id, edge_id leave_id, bool & is_swap_enter, bool & is_swap_leave);
-        bool check_well_formed();
-        void get_path(node start, node end, svector<edge_id> & path, svector<bool> & against);
         bool is_forward_edge(edge_id e_id) const;
         bool is_ancestor_of(node ancestor, node child);
+
+    public:      
+        thread_spanning_tree() {};
+        thread_spanning_tree(graph & g);
+        ~thread_spanning_tree() {};
+
+        void initialize(svector<edge_id> const & tree);
+        void get_descendants(node start, svector<node> & descendants);
+        
+        void update(edge_id enter_id, edge_id leave_id);                
+        void get_path(node start, node end, svector<edge_id> & path, svector<bool> & against);              
+        bool in_subtree_t2(node child);
+
+        bool check_well_formed();        
     };
 
 }

src/smt/spanning_tree_base.h

@@ -47,14 +47,14 @@ namespace smt {
         typedef int node;
 
     public:
-        virtual void initialize(svector<bool> const & upwards) {};
+        virtual void initialize(svector<edge_id> const & tree) {};
         virtual void get_descendants(node start, svector<node> & descendants) {};
         
-        virtual void update(edge_id enter_id, edge_id leave_id, bool & is_swap_enter, bool & is_swap_leave) {};
+        virtual void update(edge_id enter_id, edge_id leave_id) {};        
-        virtual bool check_well_formed() {UNREACHABLE(); return false;};
         virtual void get_path(node start, node end, svector<edge_id> & path, svector<bool> & against) {};
-        virtual bool is_forward_edge(edge_id e_id) const {UNREACHABLE(); return false;};
+        virtual bool in_subtree_t2(node child) {UNREACHABLE(); return false;};
-        virtual bool is_ancestor_of(node ancestor, node child) {UNREACHABLE(); return false;};
+
+        virtual bool check_well_formed() {UNREACHABLE(); return false;};
     };
 }
 

src/smt/spanning_tree_def.h

@@ -29,14 +29,15 @@ namespace smt {
     }
 
     template<typename Ext>
-    void thread_spanning_tree<Ext>::initialize(svector<bool> const & upwards) {
+    void thread_spanning_tree<Ext>::initialize(svector<edge_id> const & tree) {
+        m_tree = tree;
+
         unsigned num_nodes = m_graph.get_num_nodes();
         m_pred.resize(num_nodes);
         m_depth.resize(num_nodes);
         m_thread.resize(num_nodes);
-        m_upwards.resize(num_nodes);
         
-        node root = m_graph.get_num_nodes() - 1;
+        node root = num_nodes - 1;
         m_pred[root] = -1;
         m_depth[root] = 0;
         m_thread[root] = 0;
@@ -46,12 +47,11 @@ namespace smt {
             m_pred[i] = root;
             m_depth[i] = 1;
             m_thread[i] = i + 1;
-            m_upwards[i] = upwards[i];
         }
 
         TRACE("network_flow", {
             tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread); 
-            tout << pp_vector("Depths", m_depth) << pp_vector("Upwards", m_upwards);
+            tout << pp_vector("Depths", m_depth) << pp_vector("Tree", m_tree);
         });
     }
 
@@ -66,27 +66,18 @@ namespace smt {
         return u;
     }
 
-    template<typename Ext>
-    edge_id thread_spanning_tree<Ext>::get_edge_to_parent(node start) const {
-        SASSERT(m_pred[start] != -1);
-        edge_id id;
-        node end = m_pred[start];
-        VERIFY(m_upwards[start] ? m_graph.get_edge_id(start, end, id) : m_graph.get_edge_id(end, start, id));
-        return id;
-    }
-
     template<typename Ext>
     void thread_spanning_tree<Ext>::get_path(node start, node end, svector<edge_id> & path, svector<bool> & against) {
         node join = get_common_ancestor(start, end);
         path.reset();
         while (start != join) {
-            edge_id e_id = get_edge_to_parent(start);
+            edge_id e_id = m_tree[start];
             path.push_back(e_id);
             against.push_back(is_forward_edge(e_id));
             start = m_pred[start];
         }
         while (end != join) {
-            edge_id e_id = get_edge_to_parent(end);
+            edge_id e_id = m_tree[end];
             path.push_back(e_id);
             against.push_back(!is_forward_edge(e_id));
             end = m_pred[end];
@@ -104,11 +95,21 @@ namespace smt {
     template<typename Ext>
     void thread_spanning_tree<Ext>::get_descendants(node start, svector<node> & descendants) {
         descendants.reset();        
-        node u = start;
+        descendants.push_back(start);
-        while (m_depth[m_thread[u]] > m_depth[start]) {
+        node u = m_thread[start];
+        while (m_depth[u] > m_depth[start]) {         
             descendants.push_back(u);
             u = m_thread[u];
         }
+        
+    }
+
+    template<typename Ext>
+    bool thread_spanning_tree<Ext>::in_subtree_t2(node child) {
+        if (m_depth[child] < m_depth[m_root_t2]) {
+            return false;
+        }
+        return is_ancestor_of(m_root_t2, child);
     }
 
     template<typename Ext>
@@ -138,7 +139,7 @@ namespace smt {
                            q                         q
         */
     template<typename Ext>
-    void thread_spanning_tree<Ext>::update(edge_id enter_id, edge_id leave_id, bool & is_swap_enter, bool & is_swap_leave) {    
+    void thread_spanning_tree<Ext>::update(edge_id enter_id, edge_id leave_id) {    
         node p = m_graph.get_source(enter_id);
         node q = m_graph.get_target(enter_id);
         node u = m_graph.get_source(leave_id);
@@ -146,20 +147,14 @@ namespace smt {
         
         if (m_pred[u] == v) {
             std::swap(u, v);
-            is_swap_leave = true;
-        }
-        else {
-            is_swap_leave = false;
         }
+
         SASSERT(m_pred[v] == u);         
 
-        bool prev_upwards = false;   
         if (is_ancestor_of(v, p)) {
             std::swap(p, q);    
-            prev_upwards = true;            
         }
 
-        is_swap_enter = prev_upwards;
         SASSERT(is_ancestor_of(v, q));
         
         TRACE("network_flow", { 
@@ -167,33 +162,46 @@ namespace smt {
             tout << u << ", " << v << ") leaves\n";
         });
 
-        // Update m_pred (for nodes in the stem from q to v)
-        // Note: m_pred[v] == u
-        // Initialize m_upwards[q] = q_upwards
-
         node old_pred = m_pred[q];        
+        // Update stem nodes from q to v
         if (q != v) {
             for (node n = q; n != u; ) {
                 SASSERT(old_pred != u || n == v); // the last processed node is v
                 SASSERT(-1 != m_pred[old_pred]);
                 int next_old_pred = m_pred[old_pred];  
                 swap_order(n, old_pred);
-                std::swap(m_upwards[n], prev_upwards);
+                m_tree[old_pred] = m_tree[n];
-                prev_upwards = !prev_upwards; // flip previous version of upwards.
                 n = old_pred;
                 old_pred = next_old_pred;
             }     
         }
-        m_pred[q] = p;
+        else {
+            node x = get_final(p);
+            node y = m_thread[x];
+            node z = get_final(q);
+            node t = m_thread[get_final(v)];
+            node r = find_rev_thread(v);
+            m_thread[z] = y;
+            m_thread[x] = q;
+            m_thread[r] = t;
+        }
 
-        // m_thread were updated.
+        m_pred[q] = p;        
-        // update the depth.
+        m_tree[q] = enter_id;
+        m_root_t2 = q;
+
+        SASSERT(!in_subtree_t2(p));
+        SASSERT(in_subtree_t2(q));
+        SASSERT(!in_subtree_t2(u));
+        SASSERT(in_subtree_t2(v));
+        
+        // Update the depth.
 
         fix_depth(q, get_final(q));
 
         TRACE("network_flow", {
             tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread); 
-            tout << pp_vector("Depths", m_depth) << pp_vector("Upwards", m_upwards);
+            tout << pp_vector("Depths", m_depth) << pp_vector("Tree", m_tree);
             });
     }
 
@@ -265,6 +273,12 @@ namespace smt {
             SASSERT(roots[i] + roots.size() == 0 || roots[i] >= 0);            
         }   
 
+        for (unsigned i = 0; i < m_tree.size(); ++i) {           
+            node src = m_graph.get_source(m_tree[i]);
+            node tgt = m_graph.get_target(m_tree[i]);
+            SASSERT(m_pred[src] == tgt || m_pred[tgt] == src);            
+        }  
+
         return true;
     }