move tolerances to idol::Tolerance

hlefebvr · Aug 17, 2023 · 5f8ab4e · 5f8ab4e
1 parent 79d9c50
commit 5f8ab4e
Show file tree

Hide file tree

Showing 14 changed files with 111 additions and 35 deletions.
diff --git a/documentation/api/index.rst b/documentation/api/index.rst
@@ -10,3 +10,4 @@ API
     modeling/index
     optimizers/index
     problems/index
+    tolerances
diff --git a/documentation/api/tolerances.rst b/documentation/api/tolerances.rst
@@ -0,0 +1,4 @@
+Tolerances
+==========
+
+.. doxygennamespace:: idol::Tolerance
diff --git a/lib/include/modeling/expressions/AbstractExpr.h b/lib/include/modeling/expressions/AbstractExpr.h
@@ -115,7 +115,7 @@ class idol::impl::AbstractExpr {
     /**
      * Multiplies every term of the expression by `t_factor`.
      *
-     * If `t_factor` equals zero with tolerance `ToleranceForSparsity`, then the expression is emptied.
+     * If `t_factor` equals zero with tolerance `Tolerance::Sparsity`, then the expression is emptied.
      *
      * @param t_factor The multiplying factor.
      * @return The object itself.
@@ -222,7 +222,7 @@ template<class Key, class IteratorOutputT, class Hash, class EqualTo>
 idol::impl::AbstractExpr<Key, IteratorOutputT, Hash, EqualTo> &
 idol::impl::AbstractExpr<Key, IteratorOutputT, Hash, EqualTo>::operator*=(double t_factor) {
 
-    if (equals(t_factor, 0., ToleranceForSparsity)) {
+    if (equals(t_factor, 0., Tolerance::Sparsity)) {
         m_map.clear();
         return *this;
     }

diff --git a/lib/include/modeling/expressions/Constant.h b/lib/include/modeling/expressions/Constant.h
@@ -187,14 +187,14 @@ namespace idol {
     static std::ostream &operator<<(std::ostream &t_os, const Constant &t_coefficient) {
 
         const auto print_lin_term = [&t_os](const idol::Param &t_param, double t_coeff) {
-            if (!idol::equals(t_coeff, 1., idol::ToleranceForSparsity)) {
+            if (!idol::equals(t_coeff, 1., idol::Tolerance::Sparsity)) {
                 t_os << t_coeff << ' ';
             }
             t_os << t_param;
         };
 
         const auto print_quad_term = [&t_os](const std::pair<idol::Param, idol::Param> &t_pair, double t_coeff) {
-            if (!idol::equals(t_coeff, 1., idol::ToleranceForSparsity)) {
+            if (!idol::equals(t_coeff, 1., idol::Tolerance::Sparsity)) {
                 t_os << t_coeff << ' ';
             }
             t_os << t_pair.first << ' ' << t_pair.second;
@@ -204,7 +204,7 @@ namespace idol {
 
         bool first_term_has_been_printed = false;
 
-        if (!idol::equals(constant, 0., idol::ToleranceForSparsity)) {
+        if (!idol::equals(constant, 0., idol::Tolerance::Sparsity)) {
             t_os << constant;
             first_term_has_been_printed = true;
         }

diff --git a/lib/include/modeling/numericals.h b/lib/include/modeling/numericals.h
@@ -12,12 +12,80 @@
 namespace idol {
 
     static constexpr double Inf = 1e20;
-    static double ToleranceForSparsity = 1e-6;
-    static double ToleranceForRelativeGapMIP = 1e-4;
-    static double ToleranceForAbsoluteGapMIP = 1e-6;
-    static double ToleranceForAbsoluteGapPricing = 1e-6;
-    static double ToleranceForRelativeGapPricing = 1e-4;
-    static double ToleranceForIntegrality = 1e-5;
+
+    /**
+     * Stores the default high-level tolerances used in idol.
+     *
+     * It is possible for optimizers to have additional tolerance parameters, yet, the tolerances defined in this namespace
+     * should always be taken into account by the optimizer. Apart from the Sparsity tolerance,
+     * users can also change tolerance values at a local level (i.e., at an optimizer level) rather than at a global level.
+     */
+    namespace Tolerance {
+
+        /**
+         * **Default:** \f$ 10^{-6} \f$
+         *
+         * **Recommended range:** \f$ [ 10^{-10}, 10^{-5} ] \f$
+         *
+         * This tolerance is used when data is saved in a sparse manner.
+         * For instance, when a value close to zero should be stored
+         * or ignored when saving a primal point.
+         */
+        static double Sparsity;
+
+        /**
+         * **Default:** \f$ 10^{-4} \f$
+         *
+         * **Recommended range:** \f$ [ 0, +\infty ] \f$
+         *
+         * Used to declare optimality of a MIP solution by comparing with the current relative gap.
+         *
+         * The relative gap is computed as follows:
+         * \f[ RelativeGap := \frac{ |UB - LB| }{ 10^{-10} + |UB| }. \f]
+         */
+        static double MIPRelativeGap;
+
+        /**
+         * **Default:** \f$ 10^{-10} \f$
+         *
+         * **Recommended range:** \f$ [ 0, \infty ] \f$
+         *
+         * Used to declare optimality of a MIP solution by comparing with the current absolute gap.
+         *
+         * The absolute gap is computed as follows:
+         * \f[ AbsoluteGap := |UB - LB| \f]
+         */
+        static double MIPAbsoluteGap;
+
+        /**
+         * **Default:** \f$ 10^{-5} \f$
+         *
+         * **Recommended range:** \f$ [ 10^{-9}, 10^{-1} ] \f$
+         *
+         * Used to recognized integer values, i.e., a given value is considered integer when the closest integer point
+         * is closer than this tolerance.
+         */
+        static double Integer;
+
+        /**
+         *  **Default:** \f$ 10^{-6} \f$
+         *
+         * **Recommended range:** \f$ [ 10^{-9}, 10^{-2} ] \f$
+         *
+         * Used to characterized constraint satisfaction, i.e., a constraint is satisfied if it is not violated by a
+         * larger amount than this tolerance.
+         */
+        static double Feasibility;
+
+        /**
+         *  **Default:** \f$ 10^{-6} \f$
+         *
+         * **Recommended range:** \f$ [ 10^{-9}, 10^{-2} ] \f$
+         *
+         * Used to characterize optimality, i.e., all reduced costs must be smaller than this tolerance.
+         */
+        static double Optimality;
+    };
 
     static bool is_pos_inf(double t_value) {
         return t_value >= Inf;

diff --git a/lib/include/modeling/solutions/AbstractSolution.h b/lib/include/modeling/solutions/AbstractSolution.h
@@ -126,7 +126,7 @@ class idol::AbstractSolution {
 template<class KeyT, class CRTP>
 void idol::AbstractSolution<KeyT, CRTP>::set(const KeyT &t_key, double t_value) {
 
-    if (equals(t_value, 0., ToleranceForSparsity)) {
+    if (equals(t_value, 0., Tolerance::Sparsity)) {
         m_values.erase(t_key);
         return;
     }
@@ -209,7 +209,7 @@ CRTP &idol::AbstractSolution<KeyT, CRTP>::operator+=(const CRTP &t_rhs) {
         auto [it, success] = m_values.template emplace(key, value);
         if (!success) {
             it->second += value;
-            if (equals(it->second, 0., ToleranceForSparsity)) {
+            if (equals(it->second, 0., Tolerance::Sparsity)) {
                 m_values.erase(it);
             }
         }
@@ -220,7 +220,7 @@ CRTP &idol::AbstractSolution<KeyT, CRTP>::operator+=(const CRTP &t_rhs) {
 template<class KeyT, class CRTP>
 CRTP &idol::AbstractSolution<KeyT, CRTP>::operator*=(double t_factor) {
 
-    if (equals(t_factor, 0., ToleranceForSparsity)) {
+    if (equals(t_factor, 0., Tolerance::Sparsity)) {
         m_values.clear();
         return dynamic_cast<CRTP&>(*this);
     }

diff --git a/lib/include/optimizers/branch-and-bound/Optimizers_BranchAndBound.h b/lib/include/optimizers/branch-and-bound/Optimizers_BranchAndBound.h
@@ -47,6 +47,9 @@ class idol::Optimizers::BranchAndBound : public Algorithm {
     double m_root_node_best_bound = -Inf;
     double m_root_node_best_obj = +Inf;
 
+    double m_mip_relative_gap = Tolerance::MIPRelativeGap;
+    double m_mip_absolute_gap = Tolerance::MIPAbsoluteGap;
+
     std::optional<TreeNode> m_incumbent;
 protected:
     void build() override;
@@ -583,9 +586,9 @@ void idol::Optimizers::BranchAndBound<NodeInfoT>::log_node(LogLevel t_msg_level,
     const unsigned int id = t_node.id();
     char sign = ' ';
 
-    if (equals(objective_value, get_best_obj(), ToleranceForAbsoluteGapMIP)) {
+    if (equals(objective_value, get_best_obj(), m_mip_absolute_gap)) {
         sign = '-';
-    } else if (equals(objective_value, get_best_bound(), ToleranceForAbsoluteGapMIP)) {
+    } else if (equals(objective_value, get_best_bound(), m_mip_relative_gap)) {
         sign = '+';
     }
 

diff --git a/lib/include/optimizers/branch-and-bound/branching-rules/impls/MostInfeasbile.h b/lib/include/optimizers/branch-and-bound/branching-rules/impls/MostInfeasbile.h
@@ -80,13 +80,13 @@ double idol::BranchingRules::MostInfeasible<NodeT>::most_infeasible_score(const
     const double frac_value = fractional_part(t_node.info().primal_solution().get(t_var));
     //std::cout << t_var << " = " << t_node.info().primal_solution().get(t_var) << " has frac_value = " << std::setprecision(10) << frac_value << std::endl;
     //std::cout << "In tolerance = " << (frac_value <= ToleranceForIntegrality) << std::endl;
-    if (frac_value <= ToleranceForIntegrality) { return -Inf; }
+    if (frac_value <= Tolerance::Integer) { return -Inf; }
     return .5 - std::abs(.5 - frac_value);
 }
 
 template<class NodeT>
 bool idol::BranchingRules::MostInfeasible<NodeT>::is_integer(double t_x) {
-    return fractional_part(t_x) <= ToleranceForIntegrality;
+    return fractional_part(t_x) <= Tolerance::Integer;
 }
 
 template<class NodeT>
@@ -123,7 +123,7 @@ idol::BranchingRules::MostInfeasible<NodeT>::create_child_nodes(const Node<NodeT
 
     const auto [variable, score] = this->variable_selected_for_branching();
 
-    if (score <= ToleranceForIntegrality) {
+    if (score <= Tolerance::Integer) {
         throw Exception("Maximum infeasibility is less than ToleranceForIntegrality.");
     }
 

diff --git a/lib/include/optimizers/branch-and-bound/nodes/NodeUpdatorByBound.h b/lib/include/optimizers/branch-and-bound/nodes/NodeUpdatorByBound.h
@@ -70,7 +70,7 @@ void idol::NodeUpdatorByBound<NodeT>::update_bounds(Map<Var, double> &t_currentl
         }
 
         const double bound = t_is_lb ? m_model.get_var_lb(var) : m_model.get_var_ub(var);
-        if (!equals(result->second, bound, ToleranceForIntegrality)) {
+        if (!equals(result->second, bound, Tolerance::Integer)) {
             t_is_lb ? m_model.set_var_lb(var, result->second) : m_model.set_var_ub(var, result->second);
         }
         ++it;

diff --git a/lib/src/modeling/constraints/TempCtr.cpp b/lib/src/modeling/constraints/TempCtr.cpp
@@ -42,7 +42,7 @@ bool TempCtr::is_violated(const Solution::Primal &t_solution) const {
         case GreaterOrEqual: return lhs < rhs;
         default:;
     }
-    return equals(lhs, rhs, ToleranceForIntegrality);
+    return equals(lhs, rhs, Tolerance::Feasibility);
 }
 
 std::ostream &idol::operator<<(std::ostream& t_os, const TempCtr& t_temp_ctr) {

diff --git a/lib/src/modeling/expressions/Constant.cpp b/lib/src/modeling/expressions/Constant.cpp
@@ -9,14 +9,14 @@
 idol::Constant idol::Constant::Zero;
 
 idol::Constant::Constant(const Param &t_param, double t_value) : m_linear_terms({{t_param, t_value } }) {
-    if (equals(t_value, 0., ToleranceForSparsity)) {
+    if (equals(t_value, 0., Tolerance::Sparsity)) {
         m_linear_terms.clear();
     }
 }
 
 idol::Constant::Constant(const idol::Param &t_param_1, const idol::Param &t_param_2, double t_value)
     : m_quadratic_terms({{{t_param_1, t_param_2}, t_value}})  {
-    if (equals(t_value, 0., ToleranceForSparsity)) {
+    if (equals(t_value, 0., Tolerance::Sparsity)) {
         m_quadratic_terms.clear();
     }
 }
@@ -27,7 +27,7 @@ idol::Constant::Constant(double t_constant) : m_constant(t_constant) {
 
 void idol::Constant::set(const Param &t_param, double t_value) {
 
-    if (equals(t_value, 0., ToleranceForSparsity)) {
+    if (equals(t_value, 0., Tolerance::Sparsity)) {
         m_linear_terms.erase(t_param);
         return;
     }
@@ -50,7 +50,7 @@ double idol::Constant::get(const idol::Param &t_param_1, const idol::Param &t_pa
 
 idol::Constant &idol::Constant::operator*=(double t_factor) {
 
-    if (equals(t_factor, 0., ToleranceForSparsity)) {
+    if (equals(t_factor, 0., Tolerance::Sparsity)) {
         m_constant = 0;
         m_linear_terms.clear();
         m_quadratic_terms.clear();
@@ -111,35 +111,35 @@ idol::Constant &idol::Constant::operator-=(const Constant &t_term) {
 }
 
 void idol::Constant::insert_or_add(const Param &t_param, double t_value) {
-    if (equals(t_value, 0., ToleranceForSparsity)) {
+    if (equals(t_value, 0., Tolerance::Sparsity)) {
         return;
     }
 
     auto [it, success] = m_linear_terms.emplace(t_param, t_value);
     if (!success) {
         it->second += t_value;
-        if (equals(it->second, 0., ToleranceForSparsity)) {
+        if (equals(it->second, 0., Tolerance::Sparsity)) {
             m_linear_terms.erase(it);
         }
     }
 }
 
 void idol::Constant::insert_or_add(const idol::Param &t_param_1, const idol::Param &t_param_2, double t_value) {
-    if (equals(t_value, 0., ToleranceForSparsity)) {
+    if (equals(t_value, 0., Tolerance::Sparsity)) {
         return;
     }
 
     auto [it, success] = m_quadratic_terms.emplace(std::make_pair(t_param_1, t_param_2), t_value);
     if (!success) {
         it->second += t_value;
-        if (equals(it->second, 0., ToleranceForSparsity)) {
+        if (equals(it->second, 0., Tolerance::Sparsity)) {
             m_quadratic_terms.erase(it);
         }
     }
 }
 
 bool idol::Constant::is_zero() const {
-    return m_linear_terms.empty() && m_quadratic_terms.empty() && equals(m_constant, 0., ToleranceForSparsity);
+    return m_linear_terms.empty() && m_quadratic_terms.empty() && equals(m_constant, 0., Tolerance::Sparsity);
 }
 
 bool idol::Constant::is_numerical() const {
@@ -172,7 +172,7 @@ double idol::Constant::fix(const Solution::Dual &t_duals) const {
 
 void idol::Constant::set(const idol::Param &t_param_1, const idol::Param &t_param_2, double t_value) {
 
-    if (equals(t_value, 0., ToleranceForSparsity)) {
+    if (equals(t_value, 0., Tolerance::Sparsity)) {
         m_quadratic_terms.erase(std::make_pair( t_param_1, t_param_2 ));
         return;
     }

diff --git a/lib/src/optimizers/column-generation/Optimizers_ColumnGeneration.cpp b/lib/src/optimizers/column-generation/Optimizers_ColumnGeneration.cpp
@@ -525,8 +525,8 @@ void idol::Optimizers::ColumnGeneration::remove_artificial_variables() {
 }
 
 bool idol::Optimizers::ColumnGeneration::stopping_condition() const {
-    return get_absolute_gap() <= ToleranceForAbsoluteGapPricing
-           || get_relative_gap() <= ToleranceForRelativeGapPricing
+    return get_absolute_gap() <= Tolerance::MIPAbsoluteGap
+           || get_relative_gap() <= Tolerance::MIPRelativeGap
            || get_remaining_time() == 0;
 }
 

diff --git a/lib/src/optimizers/dantzig-wolfe/Optimizers_DantzigWolfeDecomposition.cpp b/lib/src/optimizers/dantzig-wolfe/Optimizers_DantzigWolfeDecomposition.cpp
@@ -108,7 +108,7 @@ void idol::Optimizers::DantzigWolfeDecomposition::apply_subproblem_bound_on_mast
         throw Exception("Unexpected RHS with non-numerical terms.");
     }
 
-    if (equals(t_value, current_rhs.numerical(), ToleranceForIntegrality)) { // Remove existing constraint for it is not useful anymore
+    if (equals(t_value, current_rhs.numerical(), Tolerance::Integer)) { // Remove existing constraint for it is not useful anymore
 
         subproblem.m_generation_pattern.linear().remove(it->second);
         m_master->remove(it->second);

diff --git a/lib/src/optimizers/solvers/Optimizers_GLPK.cpp b/lib/src/optimizers/solvers/Optimizers_GLPK.cpp
@@ -61,7 +61,7 @@ void idol::Optimizers::GLPK::set_var_attr(int t_index, int t_type, double t_lb,
 
     // Set bounds
     if (has_lb && has_ub) {
-        if (equals(t_lb, t_ub, ToleranceForIntegrality)) {
+        if (equals(t_lb, t_ub, Tolerance::Integer)) {
             glp_set_col_bnds(m_model, t_index, GLP_FX, t_lb, t_ub);
         } else {
             glp_set_col_bnds(m_model, t_index, GLP_DB, t_lb, t_ub);