debugging diff logic simple simplex

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

src/math/simplex/simplex_def.h

@@ -608,6 +608,7 @@ namespace simplex {
                 // optimal
                 return l_true;
             }
+            TRACE("simplex", tout << "x_i: v" << x_i << " x_j: v" << x_j << "\n";);
             var_info& vj = m_vars[x_j];
             if (x_i == null_var) {
                 if (inc && vj.m_upper_valid) {
@@ -635,8 +636,9 @@ namespace simplex {
             // 
 
             pivot(x_i, x_j, a_ij);
-            move_to_bound(x_i, inc == m.is_pos(a_ij));            
+            move_to_bound(x_i, inc != m.is_pos(a_ij));            
             SASSERT(well_formed_row(row(m_vars[x_j].m_base2row)));
+            TRACE("simplex", display(tout););
         }
         return l_true;
     }
@@ -651,22 +653,37 @@ namespace simplex {
         else {
             em.sub(vi.m_upper, vi.m_value, delta);
         }
+        TRACE("simplex", tout << "move " << (to_lower?"to_lower":"to_upper") 
+              << " v" << x << " " << em.to_string(delta) << "\n";);
         col_iterator it = M.col_begin(x), end = M.col_end(x);
         for (; it != end && is_pos(delta); ++it) {
+            //
+            // base_coeff*s + coeff*x + R = 0
+            // 
+            // to_lower    coeff > 0    base_coeff > 0     bound(s)
+            // ------------------------------------------------------
+            // T           T            T                  !to_lower
+            // T           T            F                   to_lower
+            // T           F            T                   to_lower
+            // T           F            F                  !to_lower
+            //
             var_t s = m_row2base[it.get_row().id()];
             var_info& vs = m_vars[s];
             numeral const& coeff = it.get_row_entry().m_coeff;
+            numeral const& base_coeff = vs.m_base_coeff;
             SASSERT(!m.is_zero(coeff));
-            bool inc_s = (m.is_pos(coeff) == to_lower);
+            bool base_to_lower = (m.is_pos(coeff) != m.is_pos(base_coeff)) == to_lower;
             eps_numeral const* bound = 0;
-            if (inc_s && vs.m_upper_valid) {
+            if (!base_to_lower && vs.m_upper_valid) {
                 bound = &vs.m_upper;
             }
-            else if (!inc_s && vs.m_lower_valid) {
+            else if (base_to_lower && vs.m_lower_valid) {
                 bound = &vs.m_lower;
             }
             if (bound) {
+                // |delta2*coeff| = |(bound-value)*base_coeff|
                 em.sub(*bound, vs.m_value, delta2);
+                em.mul(delta2, base_coeff, delta2);
                 em.div(delta2, coeff, delta2);
                 abs(delta2);
                 if (delta2 < delta) {
@@ -694,8 +711,11 @@ namespace simplex {
         for (; it != end; ++it) {
             var_t x = it->m_var;          
             if (x == v) continue; 
-            bool is_pos = m.is_pos(it->m_coeff);
+            bool is_pos = m.is_pos(it->m_coeff) == m.is_pos(m_vars[v].m_base_coeff);
             if ((is_pos && at_upper(x)) || (!is_pos && at_lower(x))) {
+                TRACE("simplex", tout << "v" << x << " pos: " << is_pos 
+                      << " at upper: " << at_upper(x) 
+                      << " at lower: " << at_lower(x) << "\n";);
                 continue; // variable cannot be used for improving bounds.
                 // TBD check?
             }
@@ -737,7 +757,11 @@ namespace simplex {
             var_info& vi = m_vars[s];
             numeral const& a_ij = it.get_row_entry().m_coeff;
             numeral const& a_ii = vi.m_base_coeff;
-            bool inc_s = m.is_neg(a_ij) ? inc : !inc;
+            bool inc_s = (m.is_pos(a_ii) != m.is_pos(a_ij)) ? inc : !inc;
+            TRACE("simplex", tout << "v" << x_j << " incs: " << inc_s 
+                  << " upper valid:" << vi.m_upper_valid 
+                  << " lower valid:" << vi.m_lower_valid << "\n";
+                  display_row(tout, r););
             if ((inc_s && !vi.m_upper_valid) || (!inc_s && !vi.m_lower_valid)) {
                 continue;
             }            
@@ -751,14 +775,15 @@ namespace simplex {
             if (is_neg(curr_gain)) {
                 curr_gain.neg();
             }
-            if (x_i == null_var || (curr_gain < gain) || 
+            if (x_i == null_var || (gain < curr_gain) || 
                 (is_zero(gain) && is_zero(curr_gain) && s < x_i)) {
                 x_i = s;
                 gain = curr_gain;
                 new_a_ij = a_ij;
+                TRACE("simplex", tout << "x_j v" << x_j << " x_i v" << x_i << " gain: ";
+                      tout << em.to_string(curr_gain) << "\n";);
             }
         }
-        TRACE("simplex", tout << "x_i v" << x_i << "\n";);
         return x_i;
     }
 

src/smt/theory_diff_logic.h

@@ -158,12 +158,13 @@ namespace smt {
             unsigned      m_asserted_atoms_lim;
             unsigned      m_asserted_qhead_old;
         };
+        typedef dl_graph<GExt> Graph;
 
         smt_params &                   m_params;
         arith_util                     m_util;
         arith_eq_adapter               m_arith_eq_adapter;
         theory_diff_logic_statistics   m_stats;
-        dl_graph<GExt>                 m_graph;
+        Graph                          m_graph;
         theory_var                     m_zero; // cache the variable representing the zero variable.
         int_vector                     m_scc_id;                  // Cheap equality propagation
         eq_prop_info_set               m_eq_prop_info_set;        // set of existing equality prop infos
@@ -356,9 +357,7 @@ namespace smt {
 
         void get_implied_bound_antecedents(edge_id bridge_edge, edge_id subsumed_edge, conflict_resolution & cr);
 
-        theory_var get_zero(sort* s) const { return m_zero; }
+        theory_var get_zero() const { return m_zero; }
-
-        theory_var get_zero(expr* e) const { return m_zero; }
 
         void inc_conflicts();
 

src/smt/theory_diff_logic_def.h

@@ -208,7 +208,7 @@ bool theory_diff_logic<Ext>::internalize_atom(app * n, bool gate_ctx) {
     }
     else {
         target = mk_var(lhs);
-        source = get_zero(lhs);
+        source = get_zero();
     }
     if (is_ge) {
         std::swap(target, source);
@@ -698,7 +698,7 @@ theory_var theory_diff_logic<Ext>::mk_num(app* n, rational const& r) {
     enode* e = 0;
     context& ctx = get_context();
     if (r.is_zero()) {
-        v = get_zero(n);
+        v = get_zero();
     }
     else if (ctx.e_internalized(n)) {
         e = ctx.get_enode(n);
@@ -706,7 +706,7 @@ theory_var theory_diff_logic<Ext>::mk_num(app* n, rational const& r) {
         SASSERT(v != null_theory_var);
     }
     else {
-        theory_var zero = get_zero(n);
+        theory_var zero = get_zero();
         e = ctx.mk_enode(n, false, false, true);
         v = mk_var(e);
         // internalizer is marking enodes as interpreted whenever the associated ast is a value and a constant.
@@ -1006,71 +1006,82 @@ inf_eps_rational<inf_rational> theory_diff_logic<Ext>::maximize(theory_var v) {
 
     IF_VERBOSE(1,
         for (unsigned i = 0; i < objective.size(); ++i) {
-            verbose_stream() << "Coefficient " << objective[i].second << " of theory_var " << objective[i].first << "\n";
+            verbose_stream() << "Coefficient " << objective[i].second 
+                             << " of theory_var " << objective[i].first << "\n";
         }
         verbose_stream() << "Free coefficient " << m_objective_consts[v] << "\n";);
 
-    
+    unsigned num_nodes = m_graph.get_num_nodes();
-#if 0
+    unsigned num_edges = m_graph.get_num_edges();
-    // disabled until fixed.
+    vector<dl_edge<GExt> > const& es = m_graph.get_all_edges();
-    
+    S.ensure_var(num_nodes + num_edges + m_objectives.size());
-
+    for (unsigned i = 0; i < num_nodes; ++i) {
-    // Objective coefficients now become balances
+        numeral const& a = m_graph.get_assignment(i);
-    vector<fin_numeral> balances(m_graph.get_num_nodes());
+        rational fin = a.get_rational().to_rational();
-    balances.fill(fin_numeral::zero());
+        rational inf = a.get_infinitesimal().to_rational();
-    fin_numeral sum = fin_numeral::zero();
+        mpq_inf q(fin.to_mpq(), inf.to_mpq());
-    for (unsigned i = 0; i < objective.size(); ++i) {
+        S.set_value(i, q);
-        fin_numeral balance(objective[i].second);
-        balances[objective[i].first] = balance;
-        sum += balance;
-    }    
-    // HACK: assume that v0 is always value 0
-    balances[0] = -sum;
-
-    TRACE("arith", display(tout););
-    network_flow<GExt> net_flow(m_graph, balances);
-    min_flow_result result = net_flow.min_cost();
-    SASSERT(result != UNBOUNDED);
-    if (result == OPTIMAL) {
-        numeral objective_value = net_flow.get_optimal_solution(m_objective_assignments[v], true) + numeral(m_objective_consts[v]);
-        IF_VERBOSE(1, verbose_stream() << "Optimal value of objective " << v << ": " << objective_value << std::endl;);
-        
-        DEBUG_CODE(
-        numeral initial_value = numeral(m_objective_consts[v]);
-        for (unsigned i = 0; i < objective.size(); ++i) {
-            initial_value += fin_numeral(objective[i].second) * m_graph.get_assignment(objective[i].first);
-        }
-        IF_VERBOSE(1, verbose_stream() << "Initial value of objective " << v << ": " << initial_value << std::endl;);
-        // FIXME: Network Simplex lose precisions when handling infinitesimals
-        SASSERT(objective_value >= initial_value.get_rational()););
-        vector<numeral> & current_assigments = m_objective_assignments[v];        
-        SASSERT(!current_assigments.empty());
-        // Normalize optimal assignments so that v0 is fixed to 0
-        for (unsigned i = 1; i < current_assigments.size(); ++i) {
-            current_assigments[i] -= current_assigments[0];
-        }
-
-        ast_manager& m = get_manager();
-        IF_VERBOSE(1,
-            verbose_stream() << "Optimal assigment:" << std::endl;
-            for (unsigned i = 0; i < objective.size(); ++i) {
-                theory_var v = objective[i].first;
-                verbose_stream() << mk_pp(get_enode(v)->get_owner(), m) << " |-> " << current_assigments[v] << std::endl;
-            });
-        rational r = objective_value.get_rational().to_rational();
-        rational i = objective_value.get_infinitesimal().to_rational();    
-        return inf_eps_rational<inf_rational>(inf_rational(r, i));
-    } 
-    else {
-        // Dual problem is infeasible, primal problem is unbounded
-        SASSERT(result == INFEASIBLE);
-        IF_VERBOSE(1, verbose_stream() << "Unbounded objective" << std::endl;);
-        return inf_eps_rational<inf_rational>::infinity();
     }
+    S.set_lower(get_zero(), mpq_inf(mpq(0), mpq(0)));
+    S.set_upper(get_zero(), mpq_inf(mpq(0), mpq(0)));
+    svector<unsigned> vars;
+    unsynch_mpq_manager mgr;
+    scoped_mpq_vector coeffs(mgr);
+    coeffs.push_back(mpq(1));
+    coeffs.push_back(mpq(-1));
+    coeffs.push_back(mpq(-1));
+    vars.resize(3);
+    for (unsigned i = 0; i < es.size(); ++i) {
+        dl_edge<GExt> const& e = es[i];
+        if (e.is_enabled()) {
+            unsigned base_var = num_nodes + i;
+            vars[0] = e.get_target();
+            vars[1] = e.get_source();
+            vars[2] = base_var;
+            S.add_row(base_var, 3, vars.c_ptr(), coeffs.c_ptr());
+            // t - s <= w
+            // t - s - b = 0, b >= w
+            numeral const& w = e.get_weight();
+            rational fin = w.get_rational().to_rational();
+            rational inf = w.get_infinitesimal().to_rational();
+            mpq_inf q(fin.to_mpq(),inf.to_mpq());
+            S.set_upper(base_var, q);
+        }
+    }
+    unsigned w = num_nodes + num_edges + v;
 
-#endif
+    // add objective function as row.
-    return inf_eps_rational<inf_rational>::infinity();
+    coeffs.reset();
+    vars.reset();
+    for (unsigned i = 0; i < objective.size(); ++i) {
+        coeffs.push_back(objective[i].second.to_mpq());
+        vars.push_back(objective[i].first);
+    }
+    coeffs.push_back(mpq(1));
+    vars.push_back(w);
+    S.add_row(w, vars.size(), vars.c_ptr(), coeffs.c_ptr());
 
+    TRACE("opt", S.display(tout); display(tout););
+    
+    // optimize    
+    lbool is_sat = S.make_feasible();
+    if (is_sat == l_undef) {
+        return inf_eps_rational<inf_rational>::infinity();        
+    }
+    TRACE("opt", S.display(tout); );
+    SASSERT(is_sat != l_false);
+    lbool is_fin = S.minimize(w);
+    switch (is_fin) {
+    case l_true: {
+        simplex::mpq_ext::eps_numeral const& val = S.get_value(w);
+        inf_rational r(-rational(val.first), -rational(val.second));
+        TRACE("opt", tout << r << " " << "\n"; );
+        return inf_eps_rational<inf_rational>(rational(0), r);
+    }
+    default:
+        TRACE("opt", tout << "unbounded\n"; );        
+        return inf_eps_rational<inf_rational>::infinity();        
+    }
 }
 
 template<typename Ext>