fix bugs in patching

Signed-off-by: Lev Nachmanson <levnach@hotmail.com>

src/math/lp/lar_solver.cpp

@@ -565,40 +565,6 @@ bool lar_solver::remove_from_basis(unsigned j) {
     return m_mpq_lar_core_solver.m_r_solver.remove_from_basis(j);
 }
 
-// val is the new value to be assigned to x[j]
-// return true iff can find a new basic column that would be feasible
-// after the pivoting
-bool lar_solver::remove_from_basis(unsigned basic_j, const mpq& val) {
-    SASSERT(is_base(basic_j));
-    impq del(0);
-    const auto& slv = m_mpq_lar_core_solver.m_r_solver;
-    bool grow = val < get_column_value(basic_j).x; // grow = true means that the monomial of the pivoted var has to grow
-    for (auto &c : A_r().m_rows[row_of_basic_column(basic_j)]) {
-        lpvar j = c.var();
-        if (j == basic_j) {
-            SASSERT(c.coeff().is_one());
-            continue;
-        }
-
-        bool can_pivot = column_is_free(j);
-        if (!can_pivot &&
-            ((grow && slv.monoid_can_increase(c))|| (!grow && slv.monoid_can_decrease(c)))) {            
-            if (del.is_zero())
-                del = impq(val) - get_column_value(basic_j);
-            
-            impq j_val = get_column_value(j) - del / c.coeff();
-            if (inside_bounds(j, j_val))
-                can_pivot = true;
-        }
-        
-        if (can_pivot) {
-            pivot_column_tableau(c.var(), row_of_basic_column(basic_j));
-            return true;
-        }
-    }            
-    return false;    
-}
-
 lar_term lar_solver::get_term_to_maximize(unsigned j_or_term) const {
     if (tv::is_term(j_or_term)) {
         return get_term(j_or_term);
@@ -1292,7 +1258,8 @@ void lar_solver::get_rid_of_inf_eps() {
     mpq delta = m_mpq_lar_core_solver.find_delta_for_strict_bounds(mpq(1));
     for (unsigned j = 0; j < number_of_vars(); j++) {
         auto & r = m_mpq_lar_core_solver.m_r_x[j];
-        r = impq(r.x + delta * r.y);
+        if (!r.y.is_zero())
+            r = impq(r.x + delta * r.y);
     }
 }
 

src/math/lp/lar_solver.h

@@ -260,7 +260,6 @@ class lar_solver : public column_namer {
     void update_delta_for_terms(const impq & delta, unsigned j, const vector<unsigned>&);
     void fill_vars_to_terms(vector<vector<unsigned>> & vars_to_terms);
     bool remove_from_basis(unsigned);
-    bool remove_from_basis(unsigned, const mpq&);
     lar_term get_term_to_maximize(unsigned ext_j) const;
     bool sum_first_coords(const lar_term& t, mpq & val) const;
     void collect_rounded_rows_to_fix();
@@ -362,14 +361,13 @@ public:
                       const ChangeReport& change_report) {
         if (is_base(j)) {
             TRACE("nla_solver", get_int_solver()->display_row_info(tout, row_of_basic_column(j)) << "\n";);
-            if (!remove_from_basis(j, val))
+            remove_from_basis(j);
-                return false;
         }
 
         impq ival(val);
         if (is_blocked(j, ival))
             return false;
-
+        TRACE("nla_solver", tout << "not blocked\n";);
         impq delta = get_column_value(j) - ival;
         for (const auto &c : A_r().column(j)) {
             unsigned row_index = c.var();

src/math/lp/nla_core.cpp

@@ -1215,7 +1215,7 @@ bool core::var_breaks_correct_monic_as_factor(lpvar j, const monic& m) const {
     if (!val(var(m)).is_zero())
         return true;
     
-    if (!val(j).is_zero()) // j was not zero: the new value does not matter - m must have another zero product
+    if (!val(j).is_zero()) // j was not zero: the new value does not matter - m must have another zero factor
         return false;
     // do we have another zero in m?       
     for (lpvar k : m) {
@@ -1287,29 +1287,38 @@ bool core::has_real(const monic& m) const {
             return true;
     return false;
 }
+
 // returns true if the patching is blocking
-bool core::patch_is_blocked(lpvar u, const monic& m, const lp::impq& ival) const {
+bool core::patch_is_blocked(lpvar u, const monic& patched_m, const lp::impq& ival) const {
-    SASSERT(m_to_refine.contains(m.var()));
+    if (m_cautious_patching &&
-    if (m_cautious_patching && !m_lar_solver.inside_bounds(u, ival))
+        (!m_lar_solver.inside_bounds(u, ival) || (var_is_int(u) && ival.is_int() == false))) {
+        TRACE("nla_solver", tout << "u = " << u << " blocked, for feas or integr\n";);
         return true; // blocking
+    }
     if (var_breaks_correct_monic(u)) {
         TRACE("nla_solver", tout << "u = " << u << " blocked as used in a correct monomial\n";);
         return true;
     }
     
-    bool ret = u == m.var() || (m.contains_var(u) && var_breaks_correct_monic_as_factor(u, m));
+    bool ret = u == patched_m.var() || (patched_m.contains_var(u) && var_breaks_correct_monic_as_factor(u, patched_m));
     
-    TRACE("nla_solver", tout << "u = " << u << ", m  = "; print_monic(m, tout) <<
+    TRACE("nla_solver", tout << "u = " << u << ", patched_m  = "; print_monic(patched_m, tout) <<
           "ret = " << ret << "\n";);
     
     return ret;
 }
 
-bool core::try_to_patch(lpvar k, const rational& v, const monic & m) {
+bool core::try_to_patch(lpvar patched_var, const rational& v, const monic & m) {
-    auto is_blocked = [this, k, m](lpvar u, const lp::impq& v)
+    auto is_blocked = [this, patched_var, m](lpvar u, const lp::impq& iv)        
-        { return u != k && patch_is_blocked(u, m, v); };
+        {
+            if (!m_lar_solver.inside_bounds(u, iv))
+                return true;
+            if (u == patched_var)
+                return false;
+            return patch_is_blocked(u, m, iv);
+        };
     auto change_report = [this](lpvar u) { update_to_refine_of_var(u); };
-    return m_lar_solver.try_to_patch(k, v, is_blocked, change_report);
+    return m_lar_solver.try_to_patch(patched_var, v, is_blocked, change_report);
 }
 
 bool in_power(const svector<lpvar>& vs, unsigned l) {
@@ -1348,31 +1357,36 @@ void core::patch_monomial(lpvar j) {
         return;
     }
 
+  
+    // Now we try patching the factor variables.
+    TRACE("nla_solver", tout << " trying squares\n";);
     // handle perfect squares
     if (m.vars().size() == 2 && m.vars()[0] == m.vars()[1]) {        
         rational root;
         if (v.is_perfect_square(root)) {
             lpvar k = m.vars()[0];
-            if (!var_is_int(k) && 
+            if (!var_breaks_correct_monic(k) && (try_to_patch(k, root, m) || try_to_patch(k, -root, m))) { 
-                !var_breaks_correct_monic(k) &&
+                TRACE("nla_solver", tout << "patched square\n";);
-                (try_to_patch(k, root, m) || try_to_patch(k, -root, m))
+                return;
-                ) { 
             }
         }
+        TRACE("nla_solver", tout << " cannot patch\n";);
         return;
     }
-    // We have v != abc. Let us suppose we patch b. Then b should
+
-    // be equal to v/ac = v/(abc/b) = b(v/abc)
+    // We have v != abc, but we need to have v = abc.
+    // If we patch b then b should be equal to v/ac = v/(abc/b) = b(v/abc)
     if (!v.is_zero()) {
         rational r = val(j) / v;
         SASSERT(m.is_sorted());
-    
+        TRACE("nla_solver", tout << "r = " << r << ", v = " << v << "\n";);
         for (unsigned l = 0; l < m.size(); l++) {
             lpvar k = m.vars()[l];
             if (!in_power(m.vars(), l) &&
                 !var_is_int(k) && 
                 !var_breaks_correct_monic(k) &&
                 try_to_patch(k, r * val(k), m)) { // r * val(k) gives the right value of k
+                TRACE("nla_solver", tout << "patched j " << j << "\n";);
                 SASSERT(mul_val(m) == var_val(m));
                 erase_from_to_refine(j);
                 break;
@@ -1385,13 +1399,15 @@ void core::patch_monomials_on_to_refine() {
     auto to_refine = m_to_refine.index();
     // the rest of the function might change m_to_refine, so have to copy
     unsigned sz = to_refine.size();
-    TRACE("nla_solver", tout << "sz = " << sz << "\n";);
+
     unsigned start = random();
     for (unsigned i = 0; i < sz; i++) {
         patch_monomial(to_refine[(start + i) % sz]);
         if (m_to_refine.size() == 0)
             break;
     }
+    TRACE("nla_solver", tout << "sz = " << sz << ", m_to_refine = " << m_to_refine.size() <<
+          (sz > m_to_refine.size()? " less" : "same" ) << "\n";);
 }
 
 void core::patch_monomials() {
@@ -1400,7 +1416,8 @@ void core::patch_monomials() {
     if (m_to_refine.size() == 0 || !m_nla_settings.expensive_patching()) {
         return;
     }
-    m_cautious_patching = false; //
+    NOT_IMPLEMENTED_YET();
+    m_cautious_patching = false; 
     patch_monomials_on_to_refine();
     m_lar_solver.push();
     save_tableau();