diff --git a/documentation/api/index.rst b/documentation/api/index.rst index f8d4a32d..8e071892 100644 --- 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 new file mode 100644 index 00000000..d90896c9 --- /dev/null +++ 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 index d69f48ea..3fbdf096 100644 --- 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 idol::impl::AbstractExpr & idol::impl::AbstractExpr::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 index 9df4bb2f..451bc4a0 100644 --- 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 &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 index fca21363..9203d879 100644 --- 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 index 25bf2d9b..49fea7e1 100644 --- a/lib/include/modeling/solutions/AbstractSolution.h +++ b/lib/include/modeling/solutions/AbstractSolution.h @@ -126,7 +126,7 @@ class idol::AbstractSolution { template void idol::AbstractSolution::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::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::operator+=(const CRTP &t_rhs) { template CRTP &idol::AbstractSolution::operator*=(double t_factor) { - if (equals(t_factor, 0., ToleranceForSparsity)) { + if (equals(t_factor, 0., Tolerance::Sparsity)) { m_values.clear(); return dynamic_cast(*this); } diff --git a/lib/include/optimizers/branch-and-bound/Optimizers_BranchAndBound.h b/lib/include/optimizers/branch-and-bound/Optimizers_BranchAndBound.h index 56d8e14f..3149da84 100644 --- 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 m_incumbent; protected: void build() override; @@ -583,9 +586,9 @@ void idol::Optimizers::BranchAndBound::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 index ffd7de00..d20fa4ea 100644 --- 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::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 bool idol::BranchingRules::MostInfeasible::is_integer(double t_x) { - return fractional_part(t_x) <= ToleranceForIntegrality; + return fractional_part(t_x) <= Tolerance::Integer; } template @@ -123,7 +123,7 @@ idol::BranchingRules::MostInfeasible::create_child_nodes(const Nodevariable_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 index 446f68d1..548f5f38 100644 --- 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::update_bounds(Map &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 index bde4f9eb..3f8bfe13 100644 --- 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 index 87102e56..e4cc4c32 100644 --- 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 index 7cefd592..6b59c23c 100644 --- 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 index 03158c3f..78b7bbc6 100644 --- 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 index 3c346e14..0feec8ed 100644 --- 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);