add model correction

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

src/math/simplex/model_based_opt.cpp

@@ -97,9 +97,10 @@ namespace opt {
     // e.g. t2/a2 >= t3/a3
     // then replace a3*x + t3 by t3/a3 - t2/a2 <= 0
     // 
-    bound_type model_based_opt::maximize(rational& value) {
+    inf_eps model_based_opt::maximize() {
         SASSERT(invariant());
         unsigned_vector other;
+        unsigned_vector bound_trail, bound_vars;
         while (!objective().m_vars.empty()) {
             TRACE("opt", tout << "tableau\n";);
             var v = objective().m_vars.back();
@@ -120,17 +121,79 @@ namespace opt {
 
                 mul_add(m_objective_id, - coeff/bound_coeff, bound_row_index);
                 m_rows[bound_row_index].m_alive = false;
+                bound_trail.push_back(bound_row_index);
+                bound_vars.push_back(x);
             }
             else {
-                return unbounded;
+                return inf_eps::infinity();
             }
         }
-        value = objective().m_value;
-        if (objective().m_type == t_lt) {
-            return strict;
+        //
+        // update the evaluation of variables to satisfy the bound.
+
+        update_values(bound_vars, bound_trail);
+
+        rational value = objective().m_value;
+        if (objective().m_type == t_lt) {            
+            return inf_eps(inf_rational(value, rational(-1)));
         }
         else {
-            return non_strict;
+            return inf_eps(inf_rational(value));
+        }
+    }
+
+
+    void model_based_opt::update_values(unsigned_vector const& bound_vars, unsigned_vector const& bound_trail) {
+        rational eps(0);
+        for (unsigned i = bound_trail.size(); i > 0; ) {
+            --i;
+            unsigned x = bound_vars[i];
+            row const& r = m_rows[bound_trail[i]];
+            rational val = r.m_coeff;
+            rational x_coeff;
+            vector<var> const& vars = r.m_vars;
+            for (unsigned j = 0; j < vars.size(); ++j) {
+                var const& v = vars[j];
+                if (x = v.m_id) {
+                    x_coeff = v.m_coeff;
+                }
+                else {
+                    val += m_var2value[v.m_id]*v.m_coeff;
+                }
+            }
+            SASSERT(!x_coeff.is_zero());
+            val /= -x_coeff;
+            // Adjust epsilon to be s
+            if (eps.is_zero() || (!val.is_zero() && eps > abs(val))) {
+                eps = abs(val)/rational(2);
+            }
+            if (eps.is_zero() || (!r.m_value.is_zero() && eps > abs(r.m_value))) {
+                eps = abs(r.m_value)/rational(2);
+            }
+            //
+            //
+            //     ax + t < 0
+            // <=> x < -t/a
+            // <=> x := -t/a - epsilon
+            // 
+            if (x_coeff.is_pos() && r.m_type == t_lt) {
+                val -= eps;
+            }
+            //
+            //     -ax + t < 0 
+            // <=> -ax < -t
+            // <=> -x < -t/a
+            // <=> x > t/a
+            // <=> x := t/a + epsilon
+            //
+
+            if (x_coeff.is_pos() && r.m_type == t_lt) {
+                val -= eps;
+            }
+            else if (x_coeff.is_neg() && r.m_type == t_lt) {
+                val += eps;
+            }
+            m_var2value[x] = val;
         }
     }
 
@@ -343,6 +406,10 @@ namespace opt {
         return v;
     }
 
+    rational model_based_opt::get_value(unsigned var) {
+        return m_var2value[var];
+    }
+
     void model_based_opt::set_row(unsigned row_id, vector<var> const& coeffs, rational const& c, ineq_type rel) {
         row& r = m_rows[row_id];
         rational val(c);

src/math/simplex/model_based_opt.h

@@ -23,6 +23,7 @@ Revision History:
 
 #include "util.h"
 #include "rational.h"
+#include"inf_eps_rational.h"
 
 namespace opt {
    
@@ -38,7 +39,7 @@ namespace opt {
         non_strict
     };
 
-
+    typedef inf_eps_rational<inf_rational> inf_eps;
     
     class model_based_opt {
     public:
@@ -71,8 +72,7 @@ namespace opt {
         bool invariant();
         bool invariant(unsigned index, row const& r);
 
-        row& objective() { return m_rows[0]; }
-        
+        row& objective() { return m_rows[0]; }        
 
         bool find_bound(unsigned x, unsigned& bound_index, rational& bound_coeff, unsigned_vector& other, bool is_pos);
         
@@ -83,6 +83,8 @@ namespace opt {
         void mul_add(unsigned row_id1, rational const& c, unsigned row_id2);
 
         void set_row(unsigned row_id, vector<var> const& coeffs, rational const& c, ineq_type rel);
+        
+        void update_values(unsigned_vector const& bound_vars, unsigned_vector const& bound_trail);
 
     public:
 
@@ -91,6 +93,9 @@ namespace opt {
         // add a fresh variable with value 'value'.
         unsigned add_var(rational const& value);
 
+        // retrieve updated value of variable.
+        rational get_value(unsigned var_id);
+
         // add a constraint. We assume that the constraint is 
         // satisfied under the values provided to the variables.
         void add_constraint(vector<var> const& coeffs, rational const& c, ineq_type r);
@@ -98,12 +103,13 @@ namespace opt {
         // Set the objective function (linear).
         void set_objective(vector<var> const& coeffs, rational const& c);
 
+        //
         // find a maximal value for the objective function over the current values.
         // in other words, the returned maximal value may not be globally optimal,
         // but the current evaluation of variables are used to select a local
         // optimal.
-        bound_type maximize(rational& value);
-
+        //
+        inf_eps maximize();
 
         void display(std::ostream& out) const;
         void display(std::ostream& out, row const& r) const;