diff --git a/src/math/lp/gomory.cpp b/src/math/lp/gomory.cpp index ce09bb03f4..d49c5650f2 100644 --- a/src/math/lp/gomory.cpp +++ b/src/math/lp/gomory.cpp @@ -479,6 +479,28 @@ public: return ret; } + // Returns true if the cut t >= k has efficacy at or above the configured threshold, + // i.e. the cut is worth keeping. efficacy = (k - t(x*)) / ||t||_2, where x* is the + // current LP solution. We compare efficacy^2 against threshold^2 to avoid a sqrt and + // keep the arithmetic exact until the final, deterministic double comparison. + bool gomory::cut_has_enough_efficacy(const lar_term& t, const mpq& k) { + mpq val(0); // t(x*) + mpq norm2(0); // ||t||_2^2 + for (lar_term::ival p : t) { + const mpq& a = p.coeff(); + val += a * lra.get_column_value(p.j()).x; + norm2 += a * a; + } + if (norm2.is_zero()) + return false; + mpq violation = k - val; // positive, since x* violates the cut t >= k + if (!violation.is_pos()) + return false; + double thr = lia.settings().gomory_cut_efficacy_threshold(); + // efficacy >= thr <=> violation^2 >= thr^2 * norm2 + return (violation * violation).get_double() >= thr * thr * norm2.get_double(); + } + lia_move gomory::get_gomory_cuts(unsigned num_cuts) { struct cut_result {lar_term t; mpq k; u_dependency *dep;}; vector big_cuts; @@ -519,7 +541,16 @@ public: lra.update_column_type_and_bound(j, lp::lconstraint_kind::LE, floor(lra.get_column_value(j).x), add_deps(cc.m_dep, row, j)); else if (cc.m_polarity == row_polarity::MIN) lra.update_column_type_and_bound(j, lp::lconstraint_kind::GE, ceil(lra.get_column_value(j).x), add_deps(cc.m_dep, row, j)); - + + // Option A: optionally discard cuts whose efficacy (the distance from the LP + // solution to the cut hyperplane, normalized by the coefficient norm) is too small. + // The cut is m_t >= m_k and the current LP solution violates it, so the + // violation m_k - m_t(x*) is positive. efficacy = violation / ||m_t||_2. + // To stay rational and deterministic we compare efficacy^2 against threshold^2, + // i.e. violation^2 >= threshold^2 * ||m_t||_2^2. + if (lia.settings().gomory_cut_efficacy_filter() && !cut_has_enough_efficacy(cc.m_t, cc.m_k)) + continue; + if (!is_small_cut(lia.get_term())) { big_cuts.push_back({cc.m_t, cc.m_k, cc.m_dep}); continue; diff --git a/src/math/lp/gomory.h b/src/math/lp/gomory.h index 0fdf4f8fe1..210d8bd464 100644 --- a/src/math/lp/gomory.h +++ b/src/math/lp/gomory.h @@ -29,6 +29,7 @@ namespace lp { class lar_solver& lra; unsigned_vector gomory_select_int_infeasible_vars(unsigned num_cuts); bool is_gomory_cut_target(lpvar j); + bool cut_has_enough_efficacy(const lar_term& t, const mpq& k); u_dependency* add_deps(u_dependency*, const row_strip&, lpvar); public: lia_move get_gomory_cuts(unsigned num_cuts); diff --git a/src/math/lp/lp_params_helper.pyg b/src/math/lp/lp_params_helper.pyg index 29a10c2d52..d04376f7d4 100644 --- a/src/math/lp/lp_params_helper.pyg +++ b/src/math/lp/lp_params_helper.pyg @@ -15,5 +15,7 @@ def_module_params(module_name='lp', ('lcube_flips', UINT, 16, 'maximal number of coordinate flips when repairing the rounded largest cube center, only relevant when lcube is true'), ('int_hammer_period', UINT, 4, 'period (in final_check calls) for the integer cut/cube heuristics (find_cube, hnf, gomory); a smaller value calls them more often'), ('random_hammers', BOOL, True, 'draw the periodic integer heuristic gates (find_cube, lcube, hnf, gomory, dio) at random with the same 1/period rate instead of a deterministic every-k-th-call modulus'), + ('gomory_cut_efficacy_filter', BOOL, False, 'discard a generated Gomory cut whose efficacy (distance from the LP solution to the cut hyperplane, normalized by the cut coefficient norm) is below gomory_cut_efficacy_threshold'), + ('gomory_cut_efficacy_threshold', DOUBLE, 0.01, 'minimal efficacy (normalized violation) required to keep a Gomory cut when gomory_cut_efficacy_filter is enabled'), )) diff --git a/src/math/lp/lp_settings.cpp b/src/math/lp/lp_settings.cpp index affc299788..6870299eab 100644 --- a/src/math/lp/lp_settings.cpp +++ b/src/math/lp/lp_settings.cpp @@ -48,6 +48,8 @@ void lp::lp_settings::updt_params(params_ref const& _p) { m_random_hammers = lp_p.random_hammers(); m_lcube = lp_p.lcube(); m_lcube_flips = lp_p.lcube_flips(); + m_gomory_cut_efficacy_filter = lp_p.gomory_cut_efficacy_filter(); + m_gomory_cut_efficacy_threshold = lp_p.gomory_cut_efficacy_threshold(); unsigned hammer_period = lp_p.int_hammer_period(); SASSERT(hammer_period != 0); m_int_find_cube_period = hammer_period; diff --git a/src/math/lp/lp_settings.h b/src/math/lp/lp_settings.h index bc1f2044f5..4ad4f5aba6 100644 --- a/src/math/lp/lp_settings.h +++ b/src/math/lp/lp_settings.h @@ -270,9 +270,13 @@ private: bool m_random_hammers = true; bool m_lcube = true; unsigned m_lcube_flips = 16; + bool m_gomory_cut_efficacy_filter = false; + double m_gomory_cut_efficacy_threshold = 0.01; public: bool lcube() const { return m_lcube; } unsigned lcube_flips() const { return m_lcube_flips; } + bool gomory_cut_efficacy_filter() const { return m_gomory_cut_efficacy_filter; } + double gomory_cut_efficacy_threshold() const { return m_gomory_cut_efficacy_threshold; } unsigned dio_calls_period() const { return m_dio_calls_period; } unsigned & dio_calls_period() { return m_dio_calls_period; } unsigned dio_calls_period_decrease() const { return m_dio_calls_period_decrease; }