expose models, working on network flow

Signed-off-by: Nikolaj Bjorner <nbjorner@microsoft.com>

src/api/api_opt.cpp

@@ -20,6 +20,7 @@ Revision History:
 #include"api_log_macros.h"
 #include"api_context.h"
 #include"api_util.h"
+#include"api_model.h"
 #include"opt_context.h"
 #include"cancel_eh.h"
 #include"scoped_timer.h"
@@ -122,6 +123,23 @@ extern "C" {
         Z3_CATCH_RETURN(Z3_L_UNDEF);
     }
 
+    Z3_model Z3_API Z3_optimize_get_model(Z3_context c, Z3_optimize o) {
+        Z3_TRY;
+        LOG_Z3_optimize_get_model(c, o);
+        RESET_ERROR_CODE();
+        model_ref _m;
+        to_optimize_ref(o).get_model(_m);
+        if (!_m) {
+            SET_ERROR_CODE(Z3_INVALID_USAGE);
+            RETURN_Z3(0);
+        }
+        Z3_model_ref * m_ref = alloc(Z3_model_ref); 
+        m_ref->m_model = _m;
+        mk_c(c)->save_object(m_ref);
+        RETURN_Z3(of_model(m_ref));
+        Z3_CATCH_RETURN(0);
+    }
+
     void Z3_API Z3_optimize_set_params(Z3_context c, Z3_optimize o, Z3_params p) {
         Z3_TRY;
         LOG_Z3_optimize_set_params(c, o, p);

src/api/dotnet/Optimize.cs

@@ -109,6 +109,26 @@ namespace Microsoft.Z3
             }         
         }
 
+        /// <summary>
+        /// The model of the last <c>Check</c>.
+        /// </summary>
+        /// <remarks>
+        /// The result is <c>null</c> if <c>Check</c> was not invoked before,
+        /// if its results was not <c>SATISFIABLE</c>, or if model production is not enabled.
+        /// </remarks>
+        public Model Model
+        {
+            get
+            {
+                IntPtr x = Native.Z3_optimize_get_model(Context.nCtx, NativeObject);
+                if (x == IntPtr.Zero)
+                    return null;
+                else
+                    return new Model(Context, x);
+            }
+        }
+
+	
 	public void Maximize(ArithExpr e) {
 	    Native.Z3_optimize_maximize(Context.nCtx, NativeObject, e.NativeObject);
 	}

src/api/z3_api.h

@@ -6021,6 +6021,17 @@ END_MLAPI_EXCLUDE
     Z3_lbool Z3_API Z3_optimize_check(Z3_context c, Z3_optimize o);
 
 
+    /**
+       \brief Retrieve the model for the last #Z3_optimize_check
+
+       The error handler is invoked if a model is not available because 
+       the commands above were not invoked for the given optimization 
+       solver, or if the result was \c Z3_L_FALSE.
+       
+       def_API('Z3_optimize_get_model', MODEL, (_in(CONTEXT), _in(OPTIMIZE)))
+    */
+    Z3_model Z3_API Z3_optimize_get_model(Z3_context c, Z3_optimize o);
+
     /**
        \brief Set parameters on optimization context.       
 

src/opt/opt_context.cpp

@@ -82,6 +82,10 @@ namespace opt {
         }
     }
 
+    void context::get_model(model_ref& mdl) {
+        get_solver().get_model(mdl);
+    }
+
     lbool context::execute_min_max(unsigned index, bool committed, bool is_max) {
         // HACK: reuse m_optsmt without regard for box reuse and not considering
         // use-case of lex.

src/opt/opt_context.h

@@ -75,6 +75,8 @@ namespace opt {
 
         lbool optimize();
 
+        void get_model(model_ref& m);
+
         void set_cancel(bool f);
         void reset_cancel() { set_cancel(false); }
         void cancel() { set_cancel(true); }

src/smt/network_flow.h

@@ -79,7 +79,7 @@ namespace smt {
             svector<edge_state> & m_states;
             vector<numeral> &     m_potentials;
             edge_id &             m_enter_id;
-
+            bool edge_in_tree(edge_id id) const { return m_states[id] == BASIS; }
         public: 
             pivot_rule_impl(graph & g, vector<numeral> & potentials, 
                             svector<edge_state> & states, edge_id & enter_id) 
@@ -88,46 +88,21 @@ namespace smt {
                   m_states(states),
                   m_enter_id(enter_id) {
             }
+            virtual ~pivot_rule_impl() {}
             virtual bool choose_entering_edge() = 0;
+            virtual pivot_rule rule() const = 0;
         };
         
         class first_eligible_pivot : public pivot_rule_impl {
-        private:
             edge_id m_next_edge;
-
         public:
             first_eligible_pivot(graph & g, vector<numeral> & potentials, 
                                  svector<edge_state> & states, edge_id & enter_id) : 
                 pivot_rule_impl(g, potentials, states, enter_id),
                 m_next_edge(0) {
             }
-
+            virtual bool choose_entering_edge(); 
-            bool choose_entering_edge() {
+            virtual pivot_rule rule() const { return FIRST_ELIGIBLE; }
-                TRACE("network_flow", tout << "choose_entering_edge...\n";);        
-                int num_edges = m_graph.get_num_edges();
-                for (int i = m_next_edge; i < m_next_edge + num_edges; ++i) {
-                    edge_id id = (i >= num_edges) ? (i - num_edges) : i;
-                    node src = m_graph.get_source(id);
-                    node tgt = m_graph.get_target(id);
-                    if (m_states[id] != BASIS) {
-                        numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
-                        if (cost.is_pos()) {
-                            m_enter_id = id;
-                            TRACE("network_flow", {
-                                tout << "Found entering edge " << id << " between node ";
-                                tout << src << " and node " << tgt << " with reduced cost = " << cost << "...\n";
-                            });
-                            m_next_edge = m_enter_id;
-                            if (m_next_edge >= num_edges) {
-                                m_next_edge -= num_edges;
-                            }
-                            return true;
-                        }
-                    }
-                }
-                TRACE("network_flow", tout << "Found no entering edge...\n";);
-                return false;
-            };
         };
 
         class best_eligible_pivot : public pivot_rule_impl {
@@ -136,33 +111,8 @@ namespace smt {
                                  svector<edge_state> & states, edge_id & enter_id) : 
                 pivot_rule_impl(g, potentials, states, enter_id) {
             }
-
+            virtual pivot_rule rule() const { return BEST_ELIGIBLE; }
-            bool choose_entering_edge() {
+            virtual bool choose_entering_edge();
-                TRACE("network_flow", tout << "choose_entering_edge...\n";);        
-                unsigned num_edges = m_graph.get_num_edges();
-                numeral max = numeral::zero();
-                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) {
-                        numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(i);
-                        if (cost > max) {
-                            max = cost;
-                            m_enter_id = i;
-                        }
-                    }
-                }
-                if (max.is_pos()) {
-                    TRACE("network_flow", {
-                        tout << "Found entering edge " << m_enter_id << " between node ";
-                        tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
-                        tout << " with reduced cost = " << max << "...\n";
-                    });
-                    return true;
-                }
-                TRACE("network_flow", tout << "Found no entering edge...\n";);
-                return false;
-            };
         };
 
         class candidate_list_pivot : public pivot_rule_impl {
@@ -186,93 +136,21 @@ namespace smt {
                 m_candidates(m_num_candidates) {
             }
 
-            bool choose_entering_edge() {
+            virtual pivot_rule rule() const { return CANDIDATE_LIST; }
-                if (m_current_length == 0 || m_minor_step == MINOR_STEP_LIMIT) {
-                    // Build the candidate list
-                    unsigned num_edges = m_graph.get_num_edges();
-                    numeral max = numeral::zero();
-                    m_current_length = 0;
-                    for (unsigned i = m_next_edge; i < m_next_edge + num_edges; ++i) {
-                        edge_id id = (i >= num_edges) ? i - num_edges : i;
-                        node src = m_graph.get_source(id);
-                        node tgt = m_graph.get_target(id);
-                        if (m_states[id] != BASIS) {
-                            numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
-                            if (cost.is_pos()) {
-                                m_candidates[m_current_length] = id;
-                                ++m_current_length;
-                                if (cost > max) {
-                                    max = cost;
-                                    m_enter_id = id;
-                                }
-                            }
-                            if (m_current_length >= m_num_candidates) break;
-                        }
-                    }
-                    m_next_edge = m_enter_id;
-                    m_minor_step = 1;
-                    if (max.is_pos()) {
-                        TRACE("network_flow", {
-                            tout << "Found entering edge " << m_enter_id << " between node ";
-                            tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
-                            tout << " with reduced cost = " << max << "...\n";
-                        });
-                        return true;
-                    }
-                    TRACE("network_flow", tout << "Found no entering edge...\n";);
-                    return false;
-                }
 
-                ++m_minor_step;
+            virtual bool choose_entering_edge();
-                numeral max = numeral::zero();
-                for (unsigned i = 0; i < m_current_length; ++i) {
-                    edge_id id = m_candidates[i];
-                    node src = m_graph.get_source(id);
-                    node tgt = m_graph.get_target(id);
-                    if (m_states[id] != BASIS) {
-                        numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
-                        if (cost > max) {
-                            max = cost;
-                            m_enter_id = id;
-                        }
-                        // Remove stale candidates
-                        if (!cost.is_pos()) {
-                            --m_current_length;
-                            m_candidates[i] = m_candidates[m_current_length];
-                            --i;
-                        }
-                    }
-                }
-                if (max.is_pos()) {
-                    TRACE("network_flow", {
-                        tout << "Found entering edge " << m_enter_id << " between node ";
-                        tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
-                        tout << " with reduced cost = " << max << "...\n";
-                    });
-                    return true;
-                }
-                TRACE("network_flow", tout << "Found no entering edge...\n";);
-                return false;
-            };
         };
         
         graph                m_graph;
-        spanning_tree_base * m_tree;
+        scoped_ptr<spanning_tree_base> m_tree;
-
+        scoped_ptr<pivot_rule_impl>    m_pivot;
-        // Denote supply/demand b_i on node i
+        vector<fin_numeral>  m_balances;     // Denote supply/demand b_i on node i
-        vector<fin_numeral> m_balances;
+        vector<numeral>      m_potentials;   // Duals of flows which are convenient to compute dual solutions        
-
+        vector<numeral>      m_flows;        // Basic feasible flows
-        // Duals of flows which are convenient to compute dual solutions
-        vector<numeral> m_potentials;
-
-        // Basic feasible flows
-        vector<numeral> m_flows;
-        
         svector<edge_state>  m_states;
-
         unsigned             m_step;
-
+        edge_id              m_enter_id;
-        edge_id m_enter_id, m_leave_id;
+        edge_id              m_leave_id;
         optional<numeral>    m_delta;
 
         // Initialize the network with a feasible spanning tree
@@ -290,6 +168,8 @@ namespace smt {
 
         void update_spanning_tree();
 
+        numeral get_cost() const;
+
         bool edge_in_tree(edge_id id) const;
 
         bool is_infeasible();

src/smt/network_flow_def.h

@@ -26,6 +26,128 @@ Notes:
 
 namespace smt {
 
+    template<typename Ext>    
+    bool network_flow<Ext>::first_eligible_pivot::choose_entering_edge() {
+        TRACE("network_flow", tout << "choose_entering_edge...\n";);        
+        int num_edges = m_graph.get_num_edges();
+        for (int i = m_next_edge; i < m_next_edge + num_edges; ++i) {
+            edge_id id = i % num_edges;
+            if (edge_in_tree(id)) {
+                continue;
+            }
+            node src = m_graph.get_source(id);
+            node tgt = m_graph.get_target(id);
+            numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
+            if (cost.is_pos()) {
+                m_next_edge = m_enter_id = id;
+                TRACE("network_flow", 
+                      tout << "Found entering edge " << id << " between node ";
+                      tout << src << " and node " << tgt << " with reduced cost = " << cost << "...\n";);
+                return true;
+            }
+        }
+        TRACE("network_flow", tout << "Found no entering edge...\n";);
+        return false;
+    };
+
+    template<typename Ext>    
+    bool network_flow<Ext>::best_eligible_pivot::choose_entering_edge() {
+        TRACE("network_flow", tout << "choose_entering_edge...\n";);        
+        unsigned num_edges = m_graph.get_num_edges();
+        numeral max = numeral::zero();
+        for (unsigned i = 0; i < num_edges; ++i) {
+            node src = m_graph.get_source(i);
+            node tgt = m_graph.get_target(i);
+            if (!edge_in_tree(i)) {
+                numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(i);
+                if (cost > max) {
+                    max = cost;
+                    m_enter_id = i;
+                }
+            }
+        }
+        if (max.is_pos()) {
+            TRACE("network_flow", {
+                    tout << "Found entering edge " << m_enter_id << " between node ";
+                    tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
+                    tout << " with reduced cost = " << max << "...\n";
+                });
+            return true;
+        }
+        TRACE("network_flow", tout << "Found no entering edge...\n";);
+        return false;
+    };
+
+    template<typename Ext>    
+    bool network_flow<Ext>::candidate_list_pivot::choose_entering_edge() {
+        if (m_current_length == 0 || m_minor_step == MINOR_STEP_LIMIT) {
+            // Build the candidate list
+            unsigned num_edges = m_graph.get_num_edges();
+            numeral max = numeral::zero();
+            m_current_length = 0;
+            for (unsigned i = m_next_edge; i < m_next_edge + num_edges; ++i) {
+                edge_id id = (i >= num_edges) ? i - num_edges : i;
+                node src = m_graph.get_source(id);
+                node tgt = m_graph.get_target(id);
+                if (!edge_in_tree(id)) {
+                    numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
+                    if (cost.is_pos()) {
+                        m_candidates[m_current_length] = id;
+                        ++m_current_length;
+                        if (cost > max) {
+                            max = cost;
+                            m_enter_id = id;
+                        }
+                    }
+                    if (m_current_length >= m_num_candidates) break;
+                }
+            }
+            m_next_edge = m_enter_id;
+            m_minor_step = 1;
+            if (max.is_pos()) {
+                TRACE("network_flow", {
+                        tout << "Found entering edge " << m_enter_id << " between node ";
+                        tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
+                        tout << " with reduced cost = " << max << "...\n";
+                    });
+                return true;
+            }
+            TRACE("network_flow", tout << "Found no entering edge...\n";);
+            return false;
+        }
+        
+        ++m_minor_step;
+        numeral max = numeral::zero();
+        for (unsigned i = 0; i < m_current_length; ++i) {
+            edge_id id = m_candidates[i];
+            node src = m_graph.get_source(id);
+            node tgt = m_graph.get_target(id);
+            if (!edge_in_tree(id)) {
+                numeral cost = m_potentials[src] - m_potentials[tgt] - m_graph.get_weight(id);
+                if (cost > max) {
+                    max = cost;
+                    m_enter_id = id;
+                }
+                // Remove stale candidates
+                if (!cost.is_pos()) {
+                    --m_current_length;
+                    m_candidates[i] = m_candidates[m_current_length];
+                    --i;
+                }
+            }
+        }
+        if (max.is_pos()) {
+            TRACE("network_flow", {
+                    tout << "Found entering edge " << m_enter_id << " between node ";
+                    tout << m_graph.get_source(m_enter_id) << " and node " << m_graph.get_target(m_enter_id);
+                    tout << " with reduced cost = " << max << "...\n";
+                });
+            return true;
+        }
+        TRACE("network_flow", tout << "Found no entering edge...\n";);
+        return false;
+    };
+
     template<typename Ext>
     network_flow<Ext>::network_flow(graph & g, vector<fin_numeral> const & balances) :
         m_balances(balances) {
@@ -42,7 +164,7 @@ namespace smt {
             }
         }
         TRACE("network_flow", {
-            tout << "Different logic optimization:" << std::endl;
+            tout << "Difference logic optimization:" << std::endl;
             display_dual(tout);
             tout << "Minimum cost flow:" << std::endl;
             display(tout);
@@ -52,6 +174,276 @@ namespace smt {
         m_tree = alloc(basic_spanning_tree<Ext>, m_graph);
     }
 
+
+    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_graph.get_num_nodes();
+        unsigned num_edges = m_graph.get_num_edges();
+        
+        node root = num_nodes;
+        m_graph.init_var(root);
+
+        m_potentials.resize(num_nodes + 1);     
+        m_potentials[root] = numeral::zero();
+
+        m_balances.resize(num_nodes + 1);
+        fin_numeral sum_supply = fin_numeral::zero();
+        for (unsigned i = 0; i < num_nodes; ++i) {
+            sum_supply += m_balances[i];
+        }        
+        m_balances[root] = -sum_supply;
+
+        m_flows.resize(num_nodes + num_edges);
+        m_flows.fill(numeral::zero());
+        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 tree
+        svector<edge_id> tree;
+        for (unsigned i = 0; i < num_nodes; ++i) {
+            bool is_forward = !m_balances[i].is_neg();
+            m_states[num_edges + i] = BASIS;
+            node src = is_forward ? i : root;
+            node tgt = is_forward ? root : i;            
+            m_flows[num_edges + i] = is_forward ? m_balances[i] : -m_balances[i];
+            m_potentials[i] = is_forward ? numeral::one() : -numeral::one();
+            tree.push_back(m_graph.add_edge(src, tgt, numeral::one(), explanation()));
+        }
+
+        m_tree->initialize(tree);
+
+        TRACE("network_flow", 
+              tout << pp_vector("Potentials", m_potentials);
+              tout << pp_vector("Flows", m_flows);
+              tout << "Cost: " << get_cost() << "\n";
+              tout << "Spanning tree:\n";               
+              display_spanning_tree(tout););
+        SASSERT(check_well_formed());
+    }
+
+    template<typename Ext>
+    void network_flow<Ext>::update_potentials() {        
+        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;
+        node start;
+        if (m_tree->in_subtree_t2(tgt)) {
+            change = cost;
+            start = tgt;
+        }
+        else {
+            change = -cost;
+            start = src;
+        }
+        SASSERT(m_tree->in_subtree_t2(start));
+        TRACE("network_flow", tout << "update_potentials of T_" << start << " with change = " << change << "...\n";);
+        svector<node> descendants;
+        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;
+        }
+        TRACE("network_flow", tout << pp_vector("Potentials", m_potentials););
+    }
+
+    template<typename Ext>
+    void network_flow<Ext>::update_flows() {
+        TRACE("network_flow", tout << "update_flows...\n";);
+        m_flows[m_enter_id] += *m_delta;
+        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] ? - *m_delta : *m_delta;
+        }
+        TRACE("network_flow", tout << pp_vector("Flows", m_flows););
+    }
+    
+    template<typename Ext>
+    bool network_flow<Ext>::choose_leaving_edge() {        
+        TRACE("network_flow", tout << "choose_leaving_edge...\n";);
+        node src = m_graph.get_source(m_enter_id);
+        node tgt = m_graph.get_target(m_enter_id); 
+        m_delta.set_invalid();
+        edge_id leave_id = null_edge_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];
+            if (against[i] && (!m_delta || m_flows[e_id] < *m_delta)) {
+                m_delta = m_flows[e_id];
+                leave_id = e_id;
+            }
+        }
+        m_leave_id = leave_id;
+        
+        TRACE("network_flow", 
+              tout << "Leaving edge " << m_leave_id;
+              tout << " between node " << m_graph.get_source(m_leave_id);
+              tout << " and node " << m_graph.get_target(m_leave_id) ;
+              if (m_delta) tout << " with delta = " << *m_delta;
+              tout << "\n";);            
+
+        return m_delta;
+    }
+
+    template<typename Ext>
+    void network_flow<Ext>::update_spanning_tree() {  
+        m_tree->update(m_enter_id, m_leave_id);
+    }
+
+    template<typename Ext>
+    bool network_flow<Ext>::choose_entering_edge(pivot_rule pr) {
+        if (!m_pivot || pr != m_pivot->rule()) {
+            switch (pr) {
+            case FIRST_ELIGIBLE: 
+                m_pivot = alloc(first_eligible_pivot, m_graph, m_potentials, m_states, m_enter_id);
+                break;
+            case BEST_ELIGIBLE:
+                m_pivot = alloc(best_eligible_pivot, m_graph, m_potentials, m_states, m_enter_id);
+                break;
+            case CANDIDATE_LIST:
+                m_pivot = alloc(candidate_list_pivot, m_graph, m_potentials, m_states, m_enter_id);
+                break;
+            default:
+                UNREACHABLE();
+            }
+        }
+        return m_pivot->choose_entering_edge();
+    }
+
+    // Minimize cost flows
+    template<typename Ext>
+    min_flow_result network_flow<Ext>::min_cost(pivot_rule pr) {
+        initialize();
+        while (choose_entering_edge(pr)) { 
+            bool bounded = choose_leaving_edge();
+            if (!bounded) return UNBOUNDED;
+            update_flows();
+            if (m_enter_id != m_leave_id) {
+                SASSERT(edge_in_tree(m_leave_id));
+                SASSERT(!edge_in_tree(m_enter_id));
+                m_states[m_enter_id] = BASIS;
+                m_states[m_leave_id] = LOWER;
+                update_spanning_tree();
+                update_potentials();                
+                TRACE("network_flow", 
+                      tout << "Spanning tree:\n"; 
+                      display_spanning_tree(tout);
+                      tout << "Cost: " << get_cost() << "\n";);
+                SASSERT(check_well_formed());
+            }
+        }
+        if (is_infeasible()) return INFEASIBLE;
+        TRACE("network_flow", tout << "Found optimal solution.\n";);
+        SASSERT(check_optimal());
+        return OPTIMAL;
+    }
+
+    template<typename Ext>
+    bool network_flow<Ext>::is_infeasible() {
+#if 0
+        // Flows of artificial arcs should be zero
+        unsigned num_nodes = m_graph.get_num_nodes();
+        unsigned num_edges = m_graph.get_num_edges();
+        for (unsigned i = 1; i < num_nodes; ++i) {
+            if (m_flows[num_edges - i].is_pos()) return true;
+        }
+#endif
+        return false;
+    }
+
+    // Get the optimal solution
+    template<typename Ext>
+    typename network_flow<Ext>::numeral network_flow<Ext>::get_optimal_solution(vector<numeral> & result, bool is_dual) {
+        numeral objective_value = get_cost();
+        result.reset();
+        if (is_dual) {            
+            result.append(m_potentials);                   
+        }
+        else {
+            result.append(m_flows);     
+        }
+        return objective_value;
+    }
+
+    template<typename Ext>
+    typename network_flow<Ext>::numeral network_flow<Ext>::get_cost() const {
+        numeral objective_value = numeral::zero();
+        unsigned num_edges = m_graph.get_num_edges();
+        for (unsigned i = 0; i < num_edges; ++i) {
+            if (edge_in_tree(i)) {
+                objective_value += m_flows[i].get_rational() * m_graph.get_weight(i);
+            }
+        }
+        return objective_value;
+    }
+
+    template<typename Ext>
+    bool network_flow<Ext>::edge_in_tree(edge_id id) const {
+        return m_states[id] == BASIS;
+    }
+    
+    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) {
+            SASSERT(!m_flows[i].is_neg());
+            SASSERT(edge_in_tree(i) || m_flows[i].is_zero());
+        }
+
+        unsigned num_edges = m_graph.get_num_edges();
+        for (unsigned i = 0; i < num_edges; ++i) {
+            if (edge_in_tree(i)) {
+                dl_var src = m_graph.get_source(i);
+                dl_var tgt = m_graph.get_target(i);
+                numeral weight = m_graph.get_weight(i);
+                SASSERT(m_potentials[src] - m_potentials[tgt] == weight);
+            }
+        }
+
+        return true;
+    }
+
+    template<typename Ext>
+    bool network_flow<Ext>::check_optimal() {
+        numeral total_cost = get_cost();
+        unsigned num_edges = m_graph.get_num_edges();
+
+        for (unsigned i = 0; i < num_edges; ++i) {
+            dl_var src = m_graph.get_source(i);
+            dl_var tgt = m_graph.get_target(i);
+            numeral weight = m_graph.get_weight(i);
+            SASSERT(m_potentials[src] - m_potentials[tgt] <= weight);
+        }
+
+        // m_flows are zero on non-basic edges
+        for (unsigned i = 0; i < m_flows.size(); ++i) {
+            SASSERT(edge_in_tree(i) || m_flows[i].is_zero());
+        }
+        numeral total_balance = numeral::zero();
+        for (unsigned i = 0; i < m_potentials.size(); ++i) {
+            total_balance += m_balances[i] * m_potentials[i];
+        }    
+        TRACE("network_flow", tout << "Total balance: " << total_balance << ", total cost: " << total_cost << std::endl;);
+        return total_cost == total_balance;
+    }
+
+    // display methods
+
     template<typename Ext>
     void network_flow<Ext>::display(std::ofstream & os) {
         vector<edge> const & es = m_graph.get_all_edges();
@@ -96,6 +488,7 @@ namespace smt {
         os << "(optimize)" << std::endl;
     }
 
+
     template<typename Ext>
     void network_flow<Ext>::display_dual(std::ofstream & os) {
         for (unsigned i = 0; i < m_graph.get_num_nodes(); ++i) {
@@ -115,271 +508,6 @@ namespace smt {
         os << "(optimize)" << std::endl;
     }
 
-    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_graph.get_num_nodes();
-        unsigned num_edges = m_graph.get_num_edges();
-        
-        node root = num_nodes;
-        m_graph.init_var(root);
-
-        m_potentials.resize(num_nodes + 1);     
-        m_potentials[root] = numeral::zero();
-
-        m_balances.resize(num_nodes + 1);
-        fin_numeral sum_supply = fin_numeral::zero();
-        for (unsigned i = 0; i < num_nodes; ++i) {
-            sum_supply += m_balances[i];
-        }        
-        m_balances[root] = -sum_supply;
-
-        m_flows.resize(num_nodes + num_edges);
-        m_flows.fill(numeral::zero());
-        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 tree
-        svector<edge_id> tree;
-        bool is_forward;
-        for (unsigned i = 0; i < num_nodes; ++i) {
-            is_forward = !m_balances[i].is_neg();
-            m_states[num_edges + i] = BASIS;
-            node src = is_forward ? i : root;
-            node tgt = is_forward ? root : i;            
-            m_flows[num_edges + i] = is_forward ? m_balances[i] : -m_balances[i];
-            m_potentials[i] = is_forward ? numeral::one() : -numeral::one();
-            tree.push_back(m_graph.add_edge(src, tgt, numeral::one(), explanation()));
-        }
-
-        m_tree->initialize(tree);
-
-        TRACE("network_flow", {
-                tout << pp_vector("Potentials", m_potentials, true) << pp_vector("Flows", m_flows);
-            });
-        TRACE("network_flow", tout << "Spanning tree:\n"; display_spanning_tree(tout););
-        SASSERT(check_well_formed());
-    }
-
-    template<typename Ext>
-    void network_flow<Ext>::update_potentials() {        
-        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;
-        node start;
-        if (m_tree->in_subtree_t2(tgt)) {
-            change = cost;
-            start = tgt;
-        }
-        else {
-            change = -cost;
-            start = src;
-        }
-        SASSERT(m_tree->in_subtree_t2(start));
-        TRACE("network_flow", tout << "update_potentials of T_" << start << " with change = " << change << "...\n";);
-        svector<node> descendants;
-        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;
-        }
-        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";);
-        m_flows[m_enter_id] += *m_delta;
-        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] ? - *m_delta : *m_delta;
-        }
-        TRACE("network_flow", tout << pp_vector("Flows", m_flows, true););
-    }
-    
-    template<typename Ext>
-    bool network_flow<Ext>::choose_leaving_edge() {        
-        TRACE("network_flow", tout << "choose_leaving_edge...\n";);
-        node src = m_graph.get_source(m_enter_id);
-        node tgt = m_graph.get_target(m_enter_id); 
-        m_delta.set_invalid();
-        edge_id leave_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];
-            if (against[i] && (!m_delta || m_flows[e_id] < *m_delta)) {
-                m_delta = m_flows[e_id];
-                leave_id = e_id;
-            }
-        }
-
-        if (m_delta) {
-            m_leave_id = leave_id;
-
-            TRACE("network_flow", { 
-                tout << "Found leaving edge " << m_leave_id;
-                tout << " between node " << m_graph.get_source(m_leave_id);
-                tout << " and node " << m_graph.get_target(m_leave_id) << " with delta = " << *m_delta << "...\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() {  
-        m_tree->update(m_enter_id, m_leave_id);
-    }
-
-    template<typename Ext>
-    bool network_flow<Ext>::choose_entering_edge(pivot_rule pr) {
-        pivot_rule_impl * pivot;
-        switch (pr) {
-        case FIRST_ELIGIBLE:
-            pivot = alloc(first_eligible_pivot, m_graph, m_potentials, m_states, m_enter_id);
-            break;
-        case BEST_ELIGIBLE:
-            pivot = alloc(best_eligible_pivot, m_graph, m_potentials, m_states, m_enter_id);
-            break;
-        case CANDIDATE_LIST:
-            pivot = alloc(candidate_list_pivot, m_graph, m_potentials, m_states, m_enter_id);
-            break;
-        default:
-            UNREACHABLE();
-        }
-        return pivot->choose_entering_edge();
-    }
-
-    // Minimize cost flows
-    template<typename Ext>
-    min_flow_result network_flow<Ext>::min_cost(pivot_rule pr) {
-        initialize();
-        while (choose_entering_edge(pr)) { 
-            bool bounded = choose_leaving_edge();
-            if (!bounded) return UNBOUNDED;
-            update_flows();
-            if (m_enter_id != m_leave_id) {
-				SASSERT(edge_in_tree(m_leave_id));
-				SASSERT(!edge_in_tree(m_enter_id));
-                m_states[m_enter_id] = BASIS;
-                m_states[m_leave_id] = LOWER;
-                update_spanning_tree();
-                update_potentials();                
-                TRACE("network_flow", tout << "Spanning tree:\n"; display_spanning_tree(tout););
-                SASSERT(check_well_formed());
-            }
-        }
-        if (is_infeasible()) return INFEASIBLE;
-        TRACE("network_flow", tout << "Found optimal solution.\n";);
-        SASSERT(check_optimal());
-        return OPTIMAL;
-    }
-
-    template<typename Ext>
-    bool network_flow<Ext>::is_infeasible() {
-#if 0
-        // Flows of artificial arcs should be zero
-        unsigned num_nodes = m_graph.get_num_nodes();
-        unsigned num_edges = m_graph.get_num_edges();
-        for (unsigned i = 1; i < num_nodes; ++i) {
-            if (m_flows[num_edges - i].is_pos()) return true;
-        }
-#endif
-        return false;
-    }
-
-    // Get the optimal solution
-    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();
-        unsigned num_edges = m_graph.get_num_edges();
-        for (unsigned i = 0; i < num_edges; ++i) {
-            if (m_states[i] == BASIS) 
-            {
-                objective_value += m_flows[i].get_rational() * m_graph.get_weight(i);
-            }
-        }
-        result.reset();
-        if (is_dual) {            
-            result.append(m_potentials);                   
-        }
-        else {
-            result.append(m_flows);     
-        }
-        return objective_value;
-    }
-
-    template<typename Ext>
-    bool network_flow<Ext>::edge_in_tree(edge_id id) const {
-        return m_states[id] == BASIS;
-    }
-    
-    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) {
-            SASSERT(!m_flows[i].is_neg());
-            SASSERT(m_states[i] == BASIS || m_flows[i].is_zero());
-        }
-
-        unsigned num_edges = m_graph.get_num_edges();
-        for (unsigned i = 0; i < num_edges; ++i) {
-            if (m_states[i] == BASIS) {
-                dl_var src = m_graph.get_source(i);
-                dl_var tgt = m_graph.get_target(i);
-                numeral weight = m_graph.get_weight(i);
-                SASSERT(m_potentials[src] - m_potentials[tgt] == weight);
-            }
-        }
-
-        return true;
-    }
-
-    template<typename Ext>
-    bool network_flow<Ext>::check_optimal() {
-        numeral total_cost = numeral::zero();
-        unsigned num_edges = m_graph.get_num_edges();
-        for (unsigned i = 0; i < num_edges; ++i) {
-            if (m_states[i] == BASIS) {
-                total_cost += m_flows[i].get_rational() * m_graph.get_weight(i);
-            }
-        }
-
-        for (unsigned i = 0; i < num_edges; ++i) {
-            dl_var src = m_graph.get_source(i);
-            dl_var tgt = m_graph.get_target(i);
-            numeral weight = m_graph.get_weight(i);
-            SASSERT(m_potentials[src] - m_potentials[tgt] <= weight);
-        }
-
-        // m_flows are zero on non-basic edges
-        for (unsigned i = 0; i < m_flows.size(); ++i) {
-            SASSERT(m_states[i] == BASIS || m_flows[i].is_zero());
-        }
-        numeral total_balance = numeral::zero();
-        for (unsigned i = 0; i < m_potentials.size(); ++i) {
-            total_balance += m_balances[i] * m_potentials[i];
-        }    
-        std::cout << "Total balance: " << total_balance << ", total cost: " << total_cost << std::endl;
-        return total_cost == total_balance;
-    }
-
     template<typename Ext>
     void network_flow<Ext>::display_spanning_tree(std::ofstream & os) {
         ++m_step;;
@@ -397,7 +525,7 @@ namespace smt {
         unsigned num_edges = m_graph.get_num_edges();
         for (unsigned i = 0; i < num_edges; ++i) {
             os << prefix << m_graph.get_source(i) << " -> " << prefix << m_graph.get_target(i);
-                if (m_states[i] == BASIS) {
+            if (edge_in_tree(i)) {
                 os << "[color=red,penwidth=3.0,label=\"" << m_flows[i] << "/" << m_graph.get_weight(i) << "\"];\n";
             }
             else {

src/smt/spanning_tree_base.h

@@ -25,15 +25,8 @@ Notes:
 
 namespace smt {
     template<typename TV>
-    inline std::string pp_vector(std::string const & label, TV v, bool has_header = false) {
+    inline std::string pp_vector(std::string const & label, TV v) {
         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] << " ";

src/smt/spanning_tree_def.h

@@ -50,7 +50,7 @@ namespace smt {
         }
 
         TRACE("network_flow", {
-            tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread); 
+            tout << pp_vector("Predecessors", m_pred) << pp_vector("Threads", m_thread); 
             tout << pp_vector("Depths", m_depth) << pp_vector("Tree", m_tree);
         });
     }
@@ -204,7 +204,7 @@ namespace smt {
         SASSERT(in_subtree_t2(v));
 
         TRACE("network_flow", {
-            tout << pp_vector("Predecessors", m_pred, true) << pp_vector("Threads", m_thread); 
+            tout << pp_vector("Predecessors", m_pred) << pp_vector("Threads", m_thread); 
             tout << pp_vector("Depths", m_depth) << pp_vector("Tree", m_tree);
             });
     }