From 46a97de1c894cee95a55409a9b104543e1951c13 Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Tue, 27 May 2025 13:47:22 +0200 Subject: [PATCH] Use Fractional everywhere --- ortools/glop/lp_solver.cc | 21 +++++----- ortools/glop/preprocessor.cc | 9 ++-- ortools/glop/preprocessor.h | 4 +- ortools/glop/primal_edge_norms.cc | 11 +++-- ortools/glop/reduced_costs.cc | 69 ++++++++++--------------------- ortools/glop/revised_simplex.cc | 12 +++--- ortools/glop/update_row.cc | 7 +++- ortools/glop/variable_values.cc | 2 +- ortools/lp_data/lp_parser.cc | 19 ++++++++- ortools/lp_data/lp_types.h | 6 +-- 10 files changed, 82 insertions(+), 78 deletions(-) diff --git a/ortools/glop/lp_solver.cc b/ortools/glop/lp_solver.cc index 2c9efc41c6..4e8abe93b5 100644 --- a/ortools/glop/lp_solver.cc +++ b/ortools/glop/lp_solver.cc @@ -325,7 +325,7 @@ Fractional ProblemObjectiveValue(const LinearProgram& lp, Fractional value) { // Returns the allowed error magnitude for something that should evaluate to // value under the given tolerance. Fractional AllowedError(Fractional tolerance, Fractional value) { - return tolerance * std::max(1.0, std::abs(value)); + return tolerance * std::max(Fractional(1.0), std::abs(value)); } } // namespace @@ -492,15 +492,15 @@ bool LPSolver::IsOptimalSolutionOnFacet(const LinearProgram& lp) { // primal values. // TODO(user): investigate whether to use the tolerances defined in // parameters.proto. - const double kReducedCostTolerance = 1e-9; - const double kBoundTolerance = 1e-7; + const Fractional kReducedCostTolerance = 1e-9; + const Fractional kBoundTolerance = 1e-7; const ColIndex num_cols = lp.num_variables(); for (ColIndex col(0); col < num_cols; ++col) { if (variable_statuses_[col] == VariableStatus::FIXED_VALUE) continue; const Fractional lower_bound = lp.variable_lower_bounds()[col]; const Fractional upper_bound = lp.variable_upper_bounds()[col]; const Fractional value = primal_values_[col]; - if (AreWithinAbsoluteTolerance(reduced_costs_[col], 0.0, + if (AreWithinAbsoluteTolerance(reduced_costs_[col], Fractional(0.0), kReducedCostTolerance) && (AreWithinAbsoluteTolerance(value, lower_bound, kBoundTolerance) || AreWithinAbsoluteTolerance(value, upper_bound, kBoundTolerance))) { @@ -513,7 +513,7 @@ bool LPSolver::IsOptimalSolutionOnFacet(const LinearProgram& lp) { const Fractional lower_bound = lp.constraint_lower_bounds()[row]; const Fractional upper_bound = lp.constraint_upper_bounds()[row]; const Fractional activity = constraint_activities_[row]; - if (AreWithinAbsoluteTolerance(dual_values_[row], 0.0, + if (AreWithinAbsoluteTolerance(dual_values_[row], Fractional(0.0), kReducedCostTolerance) && (AreWithinAbsoluteTolerance(activity, lower_bound, kBoundTolerance) || AreWithinAbsoluteTolerance(activity, upper_bound, kBoundTolerance))) { @@ -739,7 +739,8 @@ bool LPSolver::IsProblemSolutionConsistent( case VariableStatus::AT_UPPER_BOUND: // TODO(user): revert to an exact comparison once the bug causing this // to fail has been fixed. - if (!AreWithinAbsoluteTolerance(value, ub, 1e-7) || lb == ub) { + if (!AreWithinAbsoluteTolerance(value, ub, Fractional(1e-7)) || + lb == ub) { LogVariableStatusError(col, value, status, lb, ub); return false; } @@ -1002,7 +1003,7 @@ double LPSolver::ComputeMaxExpectedObjectiveError(const LinearProgram& lp) { double LPSolver::ComputePrimalValueInfeasibility(const LinearProgram& lp, bool* is_too_large) { - double infeasibility = 0.0; + Fractional infeasibility = 0.0; const Fractional tolerance = parameters_.solution_feasibility_tolerance(); const ColIndex num_cols = lp.num_variables(); for (ColIndex col(0); col < num_cols; ++col) { @@ -1032,7 +1033,7 @@ double LPSolver::ComputePrimalValueInfeasibility(const LinearProgram& lp, double LPSolver::ComputeActivityInfeasibility(const LinearProgram& lp, bool* is_too_large) { - double infeasibility = 0.0; + Fractional infeasibility = 0.0; int num_problematic_rows(0); const RowIndex num_rows = lp.num_constraints(); const Fractional tolerance = parameters_.solution_feasibility_tolerance(); @@ -1085,7 +1086,7 @@ double LPSolver::ComputeDualValueInfeasibility(const LinearProgram& lp, bool* is_too_large) { const Fractional allowed_error = parameters_.solution_feasibility_tolerance(); const Fractional optimization_sign = lp.IsMaximizationProblem() ? -1.0 : 1.0; - double infeasibility = 0.0; + Fractional infeasibility = 0.0; const RowIndex num_rows = lp.num_constraints(); for (RowIndex row(0); row < num_rows; ++row) { const Fractional dual_value = dual_values_[row]; @@ -1108,7 +1109,7 @@ double LPSolver::ComputeDualValueInfeasibility(const LinearProgram& lp, double LPSolver::ComputeReducedCostInfeasibility(const LinearProgram& lp, bool* is_too_large) { const Fractional optimization_sign = lp.IsMaximizationProblem() ? -1.0 : 1.0; - double infeasibility = 0.0; + Fractional infeasibility = 0.0; const ColIndex num_cols = lp.num_variables(); const Fractional tolerance = parameters_.solution_feasibility_tolerance(); for (ColIndex col(0); col < num_cols; ++col) { diff --git a/ortools/glop/preprocessor.cc b/ortools/glop/preprocessor.cc index 41e00bbcfa..b39a461f86 100644 --- a/ortools/glop/preprocessor.cc +++ b/ortools/glop/preprocessor.cc @@ -2430,7 +2430,8 @@ bool SingletonPreprocessor::IntegerSingletonColumnIsRemovable( const Fractional coefficient_ratio = coefficient / matrix_entry.coeff; // Check if coefficient_ratio is integer. if (!IsIntegerWithinTolerance( - coefficient_ratio, parameters_.solution_feasibility_tolerance())) { + coefficient_ratio, + Fractional(parameters_.solution_feasibility_tolerance()))) { return false; } } @@ -2439,7 +2440,8 @@ bool SingletonPreprocessor::IntegerSingletonColumnIsRemovable( if (IsFinite(constraint_lb)) { const Fractional lower_bound_ratio = constraint_lb / matrix_entry.coeff; if (!IsIntegerWithinTolerance( - lower_bound_ratio, parameters_.solution_feasibility_tolerance())) { + lower_bound_ratio, + Fractional(parameters_.solution_feasibility_tolerance()))) { return false; } } @@ -2448,7 +2450,8 @@ bool SingletonPreprocessor::IntegerSingletonColumnIsRemovable( if (IsFinite(constraint_ub)) { const Fractional upper_bound_ratio = constraint_ub / matrix_entry.coeff; if (!IsIntegerWithinTolerance( - upper_bound_ratio, parameters_.solution_feasibility_tolerance())) { + upper_bound_ratio, + Fractional(parameters_.solution_feasibility_tolerance()))) { return false; } } diff --git a/ortools/glop/preprocessor.h b/ortools/glop/preprocessor.h index 14c36103c0..dc171cb91d 100644 --- a/ortools/glop/preprocessor.h +++ b/ortools/glop/preprocessor.h @@ -83,13 +83,13 @@ class Preprocessor { // tolerance). bool IsSmallerWithinFeasibilityTolerance(Fractional a, Fractional b) const { return ::operations_research::IsSmallerWithinTolerance( - a, b, parameters_.solution_feasibility_tolerance()); + a, b, Fractional(parameters_.solution_feasibility_tolerance())); } bool IsSmallerWithinPreprocessorZeroTolerance(Fractional a, Fractional b) const { // TODO(user): use an absolute tolerance here to be even more defensive? return ::operations_research::IsSmallerWithinTolerance( - a, b, parameters_.preprocessor_zero_tolerance()); + a, b, Fractional(parameters_.preprocessor_zero_tolerance())); } ProblemStatus status_; diff --git a/ortools/glop/primal_edge_norms.cc b/ortools/glop/primal_edge_norms.cc index d0eeb1683e..5a39d73899 100644 --- a/ortools/glop/primal_edge_norms.cc +++ b/ortools/glop/primal_edge_norms.cc @@ -238,11 +238,14 @@ void PrimalEdgeNorms::UpdateEdgeSquaredNorms(ColIndex entering_col, const Fractional pivot = -direction[leaving_row]; DCHECK_NE(pivot, 0.0); + const ColIndex first_slack = + compact_matrix_.num_cols() - RowToColIndex(compact_matrix_.num_rows()); + // Note that this should be precise because of the call to // TestEnteringEdgeNormPrecision(). const Fractional entering_squared_norm = edge_squared_norms_[entering_col]; const Fractional leaving_squared_norm = - std::max(1.0, entering_squared_norm / Square(pivot)); + std::max(Fractional(1.0), entering_squared_norm / Square(pivot)); int stat_lower_bounded_norms = 0; const Fractional factor = 2.0 / pivot; @@ -253,7 +256,9 @@ void PrimalEdgeNorms::UpdateEdgeSquaredNorms(ColIndex entering_col, for (const ColIndex col : update_row.GetNonZeroPositions()) { const Fractional coeff = update_row.GetCoefficient(col); const Fractional scalar_product = - view.ColumnScalarProduct(col, direction_left_inverse); + col >= first_slack + ? direction_left_inverse_[col - first_slack] + : view.ColumnScalarProduct(col, direction_left_inverse); num_operations_ += view.ColumnNumEntries(col).value(); // Update the edge squared norm of this column. Note that the update @@ -288,7 +293,7 @@ void PrimalEdgeNorms::UpdateDevexWeights( const Fractional entering_norm = sqrt(PreciseSquaredNorm(direction)); const Fractional pivot_magnitude = std::abs(direction[leaving_row]); const Fractional leaving_norm = - std::max(1.0, entering_norm / pivot_magnitude); + std::max(Fractional(1.0), entering_norm / pivot_magnitude); for (const ColIndex col : update_row.GetNonZeroPositions()) { const Fractional coeff = update_row.GetCoefficient(col); const Fractional update_vector_norm = std::abs(coeff) * leaving_norm; diff --git a/ortools/glop/reduced_costs.cc b/ortools/glop/reduced_costs.cc index 4542c3ff03..9eb8bc098d 100644 --- a/ortools/glop/reduced_costs.cc +++ b/ortools/glop/reduced_costs.cc @@ -26,17 +26,12 @@ #include "ortools/glop/update_row.h" #include "ortools/glop/variables_info.h" #include "ortools/lp_data/lp_types.h" +#include "ortools/lp_data/lp_utils.h" #include "ortools/lp_data/scattered_vector.h" #include "ortools/lp_data/sparse.h" #include "ortools/util/bitset.h" #include "ortools/util/stats.h" -#ifdef OMP -#include -#endif - -#include "ortools/lp_data/lp_utils.h" - namespace operations_research { namespace glop { @@ -275,8 +270,8 @@ void ReducedCosts::PerturbCosts() { case VariableType::UPPER_AND_LOWER_BOUNDED: // Here we don't necessarily maintain the dual-feasibility of a dual // feasible solution, however we can always shift the variable to its - // other bound (because it is boxed) to restore dual-feasiblity. This is - // done by MakeBoxedVariableDualFeasible() at the end of the dual + // other bound (because it is boxed) to restore dual-feasibility. This + // is done by MakeBoxedVariableDualFeasible() at the end of the dual // phase-I algorithm. if (objective > 0.0) { cost_perturbations_[col] = magnitude; @@ -370,52 +365,30 @@ void ReducedCosts::ComputeReducedCosts() { } Fractional dual_residual_error(0.0); const ColIndex num_cols = matrix_.num_cols(); + const ColIndex first_slack = num_cols - RowToColIndex(matrix_.num_rows()); reduced_costs_.resize(num_cols, 0.0); const DenseBitRow& is_basic = variables_info_.GetIsBasicBitRow(); -#ifdef OMP - const int num_omp_threads = parameters_.num_omp_threads(); -#else - const int num_omp_threads = 1; -#endif - if (num_omp_threads == 1) { - for (ColIndex col(0); col < num_cols; ++col) { - reduced_costs_[col] = objective_[col] + cost_perturbations_[col] - - matrix_.ColumnScalarProduct( - col, basic_objective_left_inverse_.values); + for (ColIndex col(0); col < first_slack; ++col) { + reduced_costs_[col] = + objective_[col] + cost_perturbations_[col] - + matrix_.ColumnScalarProduct(col, basic_objective_left_inverse_.values); - // We also compute the dual residual error y.B - c_B. - if (is_basic.IsSet(col)) { - dual_residual_error = - std::max(dual_residual_error, std::abs(reduced_costs_[col])); - } - } - } else { -#ifdef OMP - // In the multi-threaded case, perform the same computation as in the - // single-threaded case above. - std::vector thread_local_dual_residual_error(num_omp_threads, - 0.0); - const int parallel_loop_size = num_cols.value(); -#pragma omp parallel for num_threads(num_omp_threads) - for (int i = 0; i < parallel_loop_size; i++) { - const ColIndex col(i); - reduced_costs_[col] = objective_[col] + objective_perturbation_[col] - - matrix_.ColumnScalarProduct( - col, basic_objective_left_inverse_.values); - - if (is_basic.IsSet(col)) { - thread_local_dual_residual_error[omp_get_thread_num()] = - std::max(thread_local_dual_residual_error[omp_get_thread_num()], - std::abs(reduced_costs_[col])); - } - } - // end of omp parallel for - for (int i = 0; i < num_omp_threads; i++) { + // We also compute the dual residual error y.B - c_B. + if (is_basic.IsSet(col)) { dual_residual_error = - std::max(dual_residual_error, thread_local_dual_residual_error[i]); + std::max(dual_residual_error, std::abs(reduced_costs_[col])); + } + } + for (ColIndex col(first_slack); col < num_cols; ++col) { + reduced_costs_[col] = objective_[col] + cost_perturbations_[col] - + basic_objective_left_inverse_[col - first_slack]; + + // We also compute the dual residual error y.B - c_B. + if (is_basic.IsSet(col)) { + dual_residual_error = + std::max(dual_residual_error, std::abs(reduced_costs_[col])); } -#endif // OMP } deterministic_time_ += diff --git a/ortools/glop/revised_simplex.cc b/ortools/glop/revised_simplex.cc index 92b9f9c0a0..e9b26c8029 100644 --- a/ortools/glop/revised_simplex.cc +++ b/ortools/glop/revised_simplex.cc @@ -478,7 +478,7 @@ ABSL_MUST_USE_RESULT Status RevisedSimplex::SolveInternal( double min_distance = kInfinity; const DenseRow& lower_bounds = variables_info_.GetVariableLowerBounds(); const DenseRow& upper_bounds = variables_info_.GetVariableUpperBounds(); - double cost_delta = 0.0; + Fractional cost_delta = 0.0; for (ColIndex col(0); col < num_cols_; ++col) { cost_delta += solution_primal_ray_[col] * objective_[col]; if (solution_primal_ray_[col] > 0 && upper_bounds[col] != kInfinity) { @@ -586,10 +586,12 @@ ABSL_MUST_USE_RESULT Status RevisedSimplex::SolveInternal( // infeasibility lower than its corresponding residual error. Note that // we already adapt the tolerance like this during the simplex // execution. - const Fractional primal_tolerance = std::max( - primal_residual, parameters_.primal_feasibility_tolerance()); + const Fractional primal_tolerance = + std::max(primal_residual, + Fractional(parameters_.primal_feasibility_tolerance())); const Fractional dual_tolerance = - std::max(dual_residual, parameters_.dual_feasibility_tolerance()); + std::max(dual_residual, + Fractional(parameters_.dual_feasibility_tolerance())); const Fractional primal_infeasibility = variable_values_.ComputeMaximumPrimalInfeasibility(); const Fractional dual_infeasibility = @@ -2747,7 +2749,7 @@ Status RevisedSimplex::Polish(TimeLimit* time_limit) { const auto get_diff = [this](ColIndex col, Fractional old_value, Fractional new_value) { if (col >= integrality_scale_.size() || integrality_scale_[col] == 0.0) { - return 0.0; + return Fractional(0.0); } const Fractional s = integrality_scale_[col]; return (std::abs(new_value * s - std::round(new_value * s)) - diff --git a/ortools/glop/update_row.cc b/ortools/glop/update_row.cc index cc10283d6e..84b5349911 100644 --- a/ortools/glop/update_row.cc +++ b/ortools/glop/update_row.cc @@ -304,6 +304,9 @@ void UpdateRow::ComputeUpdatesColumnWise() { non_zero_position_list_.resize(matrix_.num_cols().value()); auto* non_zeros = non_zero_position_list_.data(); + const ColIndex first_slack = + matrix_.num_cols() - RowToColIndex(matrix_.num_rows()); + const Fractional drop_tolerance = parameters_.drop_tolerance(); const auto output_coeffs = coefficient_.view(); const auto view = matrix_.view(); @@ -311,7 +314,9 @@ void UpdateRow::ComputeUpdatesColumnWise() { for (const ColIndex col : variables_info_.GetIsRelevantBitRow()) { // Coefficient of the column right inverse on the 'leaving_row'. const Fractional coeff = - view.ColumnScalarProduct(col, unit_row_left_inverse); + col >= first_slack + ? unit_row_left_inverse_[col - first_slack] + : view.ColumnScalarProduct(col, unit_row_left_inverse); // Nothing to do if 'coeff' is (almost) zero which does happen due to // sparsity. Note that it shouldn't be too bad to use a non-zero drop diff --git a/ortools/glop/variable_values.cc b/ortools/glop/variable_values.cc index 734ee5fd4f..c9ee205e19 100644 --- a/ortools/glop/variable_values.cc +++ b/ortools/glop/variable_values.cc @@ -170,7 +170,7 @@ Fractional VariableValues::ComputeSumOfPrimalInfeasibilities() const { for (ColIndex col(0); col < num_cols; ++col) { const Fractional infeasibility = GetColInfeasibility(col, values, lower_bounds, upper_bounds); - sum += std::max(0.0, infeasibility); + sum += std::max(Fractional{0.0}, infeasibility); } return sum; } diff --git a/ortools/lp_data/lp_parser.cc b/ortools/lp_data/lp_parser.cc index 2aa6879050..f41e28626a 100644 --- a/ortools/lp_data/lp_parser.cc +++ b/ortools/lp_data/lp_parser.cc @@ -235,6 +235,21 @@ bool LPParser::ParseConstraint(StringPiece constraint) { return true; } +namespace { + +template +bool SimpleAtoFractional(absl::string_view str, T* value) { + if constexpr (std::is_same_v) { + return absl::SimpleAtod(str, value); + } else if constexpr (std::is_same_v) { + return absl::SimpleAtof(str, value); + } else { + static_assert(false, "Unsupported fractional type"); + return false; + } +} +} // namespace + bool LPParser::SetVariableBounds(ColIndex col, Fractional lb, Fractional ub) { if (bounded_variables_.find(col) == bounded_variables_.end()) { // The variable was not bounded yet, thus reset its bounds. @@ -250,7 +265,7 @@ bool LPParser::SetVariableBounds(ColIndex col, Fractional lb, Fractional ub) { } TokenType ConsumeToken(StringPiece* sp, std::string* consumed_name, - double* consumed_coeff) { + Fractional* consumed_coeff) { DCHECK(consumed_name != nullptr); DCHECK(consumed_coeff != nullptr); // We use LazyRE2 everywhere so that all the patterns are just compiled once @@ -305,7 +320,7 @@ TokenType ConsumeToken(StringPiece* sp, std::string* consumed_name, static const LazyRE2 kValuePattern = { R"(\s*([0-9]*\.?[0-9]+([eE][-+]?[0-9]+)?))"}; if (RE2::Consume(sp, *kValuePattern, &coeff)) { - if (!absl::SimpleAtod(coeff, consumed_coeff)) { + if (!SimpleAtoFractional(coeff, consumed_coeff)) { // Note: If absl::SimpleAtod(), Consume(), and kValuePattern are correct, // this should never happen. LOG(ERROR) << "Text: " << coeff << " was matched by RE2 to be " diff --git a/ortools/lp_data/lp_types.h b/ortools/lp_data/lp_types.h index 1e45d59301..837a0a8ec9 100644 --- a/ortools/lp_data/lp_types.h +++ b/ortools/lp_data/lp_types.h @@ -81,13 +81,13 @@ static inline double ToDouble(double f) { return f; } typedef double Fractional; // Range max for type Fractional. DBL_MAX for double for example. -constexpr double kRangeMax = std::numeric_limits::max(); +constexpr Fractional kRangeMax = std::numeric_limits::max(); // Infinity for type Fractional. -constexpr double kInfinity = std::numeric_limits::infinity(); +constexpr Fractional kInfinity = std::numeric_limits::infinity(); // Epsilon for type Fractional, i.e. the smallest e such that 1.0 + e != 1.0 . -constexpr double kEpsilon = std::numeric_limits::epsilon(); +constexpr Fractional kEpsilon = std::numeric_limits::epsilon(); // Returns true if the given value is finite, that means for a double: // not a NaN and not +/- infinity.