add option to control epsilon #7791

#7791 reports on using model values during lex optimization that break soft constraints.
This is an artifact of using optimization where optimal values can be arbitrarily close to a rational.
In a way it is by design, but we give the user now an option to control the starting point for epsilon when converting infinitesimals into rationals.

src/math/lp/lar_solver.cpp

@@ -1523,7 +1523,7 @@ namespace lp {
         if (!m_imp->m_columns_with_changed_bounds.empty())
             return false;
 
-        m_imp->m_delta = get_core_solver().find_delta_for_strict_bounds(mpq(1));
+        m_imp->m_delta = get_core_solver().find_delta_for_strict_bounds(m_imp->m_settings.m_epsilon);
         unsigned j;
         unsigned n = get_core_solver().r_x().size();
         do {
@@ -1545,7 +1545,7 @@ namespace lp {
     }
 
     void lar_solver::get_model_do_not_care_about_diff_vars(std::unordered_map<lpvar, mpq>& variable_values) const {
-        mpq delta = get_core_solver().find_delta_for_strict_bounds(mpq(1));
+        mpq delta = get_core_solver().find_delta_for_strict_bounds(m_imp->m_settings.m_epsilon);
         for (unsigned i = 0; i < get_core_solver().r_x().size(); i++) {
             const impq& rp = get_core_solver().r_x(i);
             variable_values[i] = rp.x + delta * rp.y;
@@ -1569,7 +1569,7 @@ namespace lp {
         }
         if (y_is_zero)
             return;
-        mpq delta = get_core_solver().find_delta_for_strict_bounds(mpq(1));
+        mpq delta = get_core_solver().find_delta_for_strict_bounds(m_imp->m_settings.m_epsilon);
         for (unsigned j = 0; j < number_of_vars(); j++) {
             auto& v = get_core_solver().r_x(j);
             if (!v.y.is_zero()) {

src/math/lp/lp_settings.cpp

@@ -34,6 +34,8 @@ void lp::lp_settings::updt_params(params_ref const& _p) {
     report_frequency = p.arith_rep_freq();
     m_simplex_strategy = static_cast<lp::simplex_strategy_enum>(p.arith_simplex_strategy());
     m_nlsat_delay = p.arith_nl_delay();
+    auto eps = p.arith_epsilon();
+    m_epsilon = rational(std::max(1, (int)(100000*eps)), 100000); 
     m_dio = lp_p.dio();
     m_dio_enable_gomory_cuts = lp_p.dio_cuts_enable_gomory();
     m_dio_enable_hnf_cuts = lp_p.dio_cuts_enable_hnf();

src/math/lp/lp_settings.h

@@ -242,6 +242,7 @@ private:
 public:
     unsigned         limit_on_rows_for_hnf_cutter = 75;
     unsigned         limit_on_columns_for_hnf_cutter = 150;
+    mpq              m_epsilon = mpq(1);
 private:
     unsigned         m_nlsat_delay = 0;
     bool             m_enable_hnf = true;

src/params/smt_params_helper.pyg

@@ -86,6 +86,7 @@ def_module_params(module_name='smt',
 			  ('arith.nl.cross_nested', BOOL, True, 'enable cross-nested consistency checking'),
 			  ('arith.nl.log', BOOL, False, 'Log lemmas sent to nra solver'),
                           ('arith.propagate_eqs', BOOL, True, 'propagate (cheap) equalities'),
+			  ('arith.epsilon', DOUBLE, 1.0, 'initial value of epsilon used for model generation of infinitesimals'),
                           ('arith.propagation_mode', UINT, 1, '0 - no propagation, 1 - propagate existing literals, 2 - refine finite bounds'),
                           ('arith.branch_cut_ratio', UINT, 2, 'branch/cut ratio for linear integer arithmetic'),
                           ('arith.int_eq_branch', BOOL, False, 'branching using derived integer equations'),

src/params/theory_arith_params.cpp

@@ -41,6 +41,8 @@ void theory_arith_params::updt_params(params_ref const & _p) {
     m_nl_arith_propagate_linear_monomials = p.arith_nl_propagate_linear_monomials();
     m_nl_arith_optimize_bounds = p.arith_nl_optimize_bounds();
     m_nl_arith_cross_nested = p.arith_nl_cross_nested();
+    auto eps = p.arith_epsilon();
+    m_arith_epsilon = rational(std::max(1, (int)(100000*eps)), 100000); 
 
     arith_rewriter_params ap(_p);
     m_arith_eq2ineq = ap.eq2ineq();
@@ -99,4 +101,5 @@ void theory_arith_params::display(std::ostream & out) const {
     DISPLAY_PARAM(m_nl_arith_cross_nested);
     DISPLAY_PARAM(m_arith_validate);
     DISPLAY_PARAM(m_arith_dump_lemmas);
+    DISPLAY_PARAM(m_arith_epsilon);
 }

src/params/theory_arith_params.h

@@ -20,6 +20,7 @@ Revision History:
 
 #include<climits>
 #include "util/params.h"
+#include "util/rational.h"
 
 enum class arith_solver_id {
     AS_NO_ARITH,              // 0
@@ -76,6 +77,7 @@ struct theory_arith_params {
     unsigned                m_arith_branch_cut_ratio = 2;
     bool                    m_arith_int_eq_branching = false;
     bool                    m_arith_enum_const_mod = false;
+    rational                m_arith_epsilon = rational::one();
 
     bool                    m_arith_gcd_test = true;
     bool                    m_arith_eager_gcd = false;

src/smt/theory_arith_core.h

@@ -3174,7 +3174,8 @@ namespace smt {
 
     template<typename Ext>
     void theory_arith<Ext>::compute_epsilon() {
-        m_epsilon = numeral(1);
+        auto eps = ctx.get_fparams().m_arith_epsilon;        
+        m_epsilon = numeral(eps);
         theory_var num = get_num_vars();
         for (theory_var v = 0; v < num; v++) {
             bound * l = lower(v);