always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Update pdlp

Browse Source
This commit is contained in:
Corentin Le Molgat
2022-03-02 22:09:57 +01:00
parent 35ef175cac
commit 1eba970bf3
14 changed files with 229 additions and 102 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
ortools/linear_solver/linear_solver.cc
View File

@@ -369,8 +369,8 @@ extern MPSolverInterface* BuildGLPKInterface(bool mip, MPSolver* const solver);
#endif
extern MPSolverInterface* BuildBopInterface(MPSolver* const solver);
extern MPSolverInterface* BuildGLOPInterface(MPSolver* const solver);
extern MPSolverInterface* BuildSatInterface(MPSolver* const solver);
extern MPSolverInterface* BuildPdlpInterface(MPSolver* const solver);
extern MPSolverInterface* BuildSatInterface(MPSolver* const solver);
#if defined(USE_SCIP)
extern MPSolverInterface* BuildSCIPInterface(MPSolver* const solver);
#endif
@@ -378,8 +378,6 @@ extern MPSolverInterface* BuildGurobiInterface(bool mip,
MPSolver* const solver);
#if defined(USE_CPLEX)
extern MPSolverInterface* BuildCplexInterface(bool mip, MPSolver* const solver);
extern MPSolverInterface* BuildGLOPInterface(MPSolver* const solver);
#endif
#if defined(USE_XPRESS)
extern MPSolverInterface* BuildXpressInterface(bool mip,
@@ -392,12 +390,12 @@ MPSolverInterface* BuildSolverInterface(MPSolver* const solver) {
switch (solver->ProblemType()) {
case MPSolver::BOP_INTEGER_PROGRAMMING:
return BuildBopInterface(solver);
case MPSolver::SAT_INTEGER_PROGRAMMING:
return BuildSatInterface(solver);
case MPSolver::GLOP_LINEAR_PROGRAMMING:
return BuildGLOPInterface(solver);
case MPSolver::PDLP_LINEAR_PROGRAMMING:
return BuildPdlpInterface(solver);
case MPSolver::SAT_INTEGER_PROGRAMMING:
return BuildSatInterface(solver);
#if defined(USE_GLPK)
case MPSolver::GLPK_LINEAR_PROGRAMMING:
return BuildGLPKInterface(false, solver);

3
ortools/linear_solver/linear_solver.h
View File

@@ -569,7 +569,8 @@ class MPSolver {
return solver == MPModelRequest::GLOP_LINEAR_PROGRAMMING ||
solver == MPModelRequest::GUROBI_LINEAR_PROGRAMMING ||
solver == MPModelRequest::GUROBI_MIXED_INTEGER_PROGRAMMING ||
solver == MPModelRequest::SAT_INTEGER_PROGRAMMING;
solver == MPModelRequest::SAT_INTEGER_PROGRAMMING ||
solver == MPModelRequest::PDLP_LINEAR_PROGRAMMING;
}
/// Exports model to protocol buffer.

14
ortools/linear_solver/pdlp_interface.cc
View File

@@ -107,6 +107,12 @@ MPSolver::ResultStatus PdlpInterface::Solve(const MPSolverParameters& param) {
ExtractModel();
SetParameters(param);
if (quiet_) {
parameters_.set_verbosity_level(0);
} else {
parameters_.set_verbosity_level(3);
}
solver_->SetSolverSpecificParametersAsString(
solver_->solver_specific_parameter_string_);
@@ -117,8 +123,6 @@ MPSolver::ResultStatus PdlpInterface::Solve(const MPSolverParameters& param) {
static_cast<double>(solver_->time_limit()) / 1000.0);
}
// TODO(user): Adjust logging level based on quiet_.
// Mark variables and constraints as extracted.
for (int i = 0; i < solver_->variables_.size(); ++i) {
set_variable_as_extracted(i, true);
@@ -166,12 +170,8 @@ MPSolver::ResultStatus PdlpInterface::Solve(const MPSolverParameters& param) {
absl::optional<MPSolutionResponse> PdlpInterface::DirectlySolveProto(
const MPModelRequest& request, std::atomic<bool>* interrupt) {
if (interrupt != nullptr) {
return absl::nullopt;
}
absl::StatusOr<MPSolutionResponse> response =
PdlpSolveProto(request, /*relax_integer_variables=*/true);
PdlpSolveProto(request, /*relax_integer_variables=*/true, interrupt);
if (!response.ok()) {
LOG(ERROR) << "Unexpected error solving with PDLP: " << response.status();

19
ortools/linear_solver/pdlp_proto_solver.cc
View File

@@ -13,6 +13,7 @@
#include "ortools/linear_solver/pdlp_proto_solver.h"
#include <atomic>
#include <optional>
#include <string>
#include <utility>
@@ -33,8 +34,14 @@
namespace operations_research {
absl::StatusOr<MPSolutionResponse> PdlpSolveProto(
const MPModelRequest& request, const bool relax_integer_variables) {
const MPModelRequest& request, const bool relax_integer_variables,
const std::atomic<bool>* interrupt_solve) {
pdlp::PrimalDualHybridGradientParams params;
if (request.enable_internal_solver_output()) {
params.set_verbosity_level(3);
} else {
params.set_verbosity_level(0);
}
MPSolutionResponse error_response;
if (!ProtobufTextFormatMergeFromString(request.solver_specific_parameters(),
@@ -43,11 +50,14 @@ absl::StatusOr<MPSolutionResponse> PdlpSolveProto(
MPSolverResponseStatus::MPSOLVER_MODEL_INVALID_SOLVER_PARAMETERS);
return error_response;
}
if (interrupt_solve != nullptr && interrupt_solve->load() == true) {
error_response.set_status(MPSolverResponseStatus::MPSOLVER_NOT_SOLVED);
return error_response;
}
if (request.has_solver_time_limit_seconds()) {
params.mutable_termination_criteria()->set_time_sec_limit(
request.solver_time_limit_seconds());
}
// TODO(user): Adjust logging level based on enable_internal_solver_output.
const absl::optional<LazyMutableCopy<MPModelProto>> optional_model =
ExtractValidMPModelOrPopulateResponseStatus(request, &error_response);
@@ -66,7 +76,7 @@ absl::StatusOr<MPSolutionResponse> PdlpSolveProto(
const double objective_scaling_factor = qp.objective_scaling_factor;
pdlp::SolverResult pdhg_result =
pdlp::PrimalDualHybridGradient(std::move(qp), params);
pdlp::PrimalDualHybridGradient(std::move(qp), params, interrupt_solve);
// PDLP's statuses don't map very cleanly to MPSolver statuses. Do the best
// we can for now.
@@ -81,6 +91,9 @@ absl::StatusOr<MPSolutionResponse> PdlpSolveProto(
case pdlp::TERMINATION_REASON_PRIMAL_INFEASIBLE:
response.set_status(MPSOLVER_INFEASIBLE);
break;
case pdlp::TERMINATION_REASON_INTERRUPTED_BY_USER:
response.set_status(MPSOLVER_CANCELLED_BY_USER);
break;
default:
response.set_status(MPSOLVER_NOT_SOLVED);
}

27
ortools/linear_solver/pdlp_proto_solver.h
View File

@@ -14,22 +14,31 @@
#ifndef OR_TOOLS_LINEAR_SOLVER_PDLP_PROTO_SOLVER_H_
#define OR_TOOLS_LINEAR_SOLVER_PDLP_PROTO_SOLVER_H_
#include <atomic>
#include "absl/status/statusor.h"
#include "ortools/linear_solver/linear_solver.pb.h"
namespace operations_research {
// Uses pdlp::PrimalDualHybridGradient to solve the problem specified by the
// MPModelRequest. If relax_integer_variables is true, integrality constraints
// are relaxed before solving. If false, integrality constraints result in an
// error. The solver_specific_info field in the MPSolutionResponse contains a
// serialized SolveLog. Users of this interface should be aware of the size
// limitations of MPModelProto (see, e.g., large_linear_program.proto). Returns
// an error if the conversion from MPModelProto to pdlp::QuadraticProgram fails.
// The lack of an error does not imply success. Check the SolveLog's
// termination_reason for more refined status details.
// MPModelRequest. Users of this interface should be aware of the size
// limitations of MPModelProto (see, e.g., large_linear_program.proto).
//
// The optional interrupt_solve can be used to interrupt the solve early. The
// solver will periodically check its value and stop if it holds true.
//
// If relax_integer_variables is true, integrality constraints are relaxed
// before solving. If false, integrality constraints result in an error. The
// solver_specific_info field in the MPSolutionResponse contains a serialized
// SolveLog.
//
// Returns an error if the conversion from MPModelProto to
// pdlp::QuadraticProgram fails. The lack of an error does not imply success.
// Check the SolveLog's termination_reason for more refined status details.
absl::StatusOr<MPSolutionResponse> PdlpSolveProto(
const MPModelRequest& request, bool relax_integer_variables = false);
const MPModelRequest& request, bool relax_integer_variables = false,
const std::atomic<bool>* interrupt_solve = nullptr);
} // namespace operations_research

4
ortools/linear_solver/python/linear_solver.i
View File

@@ -263,13 +263,13 @@ PY_CONVERT(MPVariable);
// Expose the MPSolver::OptimizationProblemType enum.
%unignore operations_research::MPSolver::OptimizationProblemType;
%unignore operations_research::MPSolver::GLOP_LINEAR_PROGRAMMING;
%unignore operations_research::MPSolver::CLP_LINEAR_PROGRAMMING;
%unignore operations_research::MPSolver::GLOP_LINEAR_PROGRAMMING;
%unignore operations_research::MPSolver::GLPK_LINEAR_PROGRAMMING;
%unignore operations_research::MPSolver::PDLP_LINEAR_PROGRAMMING;
%unignore operations_research::MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING;
%unignore operations_research::MPSolver::CBC_MIXED_INTEGER_PROGRAMMING;
%unignore operations_research::MPSolver::GLPK_MIXED_INTEGER_PROGRAMMING;
%unignore operations_research::MPSolver::SCIP_MIXED_INTEGER_PROGRAMMING;
%unignore operations_research::MPSolver::BOP_INTEGER_PROGRAMMING;
%unignore operations_research::MPSolver::SAT_INTEGER_PROGRAMMING;
// These aren't unit tested, as they only run on machines with a Gurobi license.

142
ortools/pdlp/primal_dual_hybrid_gradient.cc
View File

@@ -237,6 +237,7 @@ class Solver {
// Zero is used if initial_solution is nullopt.
SolverResult Solve(absl::optional<PrimalAndDualSolution> initial_solution,
const std::atomic<bool>* interrupt_solve,
IterationStatsCallback iteration_stats_callback);
private:
@@ -809,8 +810,9 @@ double Solver::ComputeNewPrimalWeight() const {
const double new_primal_weight =
std::exp(smoothing_param * std::log(unsmoothed_new_primal_weight) +
(1.0 - smoothing_param) * std::log(primal_weight_));
VLOG(4) << "New computed primal weight is " << new_primal_weight
<< " at iteration " << iterations_completed_;
LOG_IF(INFO, params_.verbosity_level() >= 4)
<< "New computed primal weight is " << new_primal_weight
<< " at iteration " << iterations_completed_;
return new_primal_weight;
}
@@ -879,13 +881,15 @@ SolverResult Solver::ConstructSolverResult(VectorXd primal_solution,
.bound_norms = original_bound_norms_});
}
VLOG(1) << "Termination reason: "
<< TerminationReason_Name(termination_reason);
VLOG(1) << "Solution point type: " << PointType_Name(output_type);
VLOG(1) << "Final solution stats:";
VLOG(1) << IterationStatsLabelString();
VLOG(1) << ToString(stats, params_.termination_criteria(),
original_bound_norms_, solve_log.solution_type());
if (params_.verbosity_level() >= 1) {
LOG(INFO) << "Termination reason: "
<< TerminationReason_Name(termination_reason);
LOG(INFO) << "Solution point type: " << PointType_Name(output_type);
LOG(INFO) << "Final solution stats:";
LOG(INFO) << IterationStatsLabelString();
LOG(INFO) << ToString(stats, params_.termination_criteria(),
original_bound_norms_, solve_log.solution_type());
}
solve_log.set_iteration_count(stats.iteration_number());
solve_log.set_termination_reason(termination_reason);
solve_log.set_solve_time_sec(stats.cumulative_time_sec());
@@ -1071,7 +1075,7 @@ Solver::UpdateIterationStatsAndCheckTermination(
}
constexpr int kLogEvery = 15;
static std::atomic_int log_counter{0};
if (VLOG_IS_ON(4)) {
if (params_.verbosity_level() >= 4) {
if (log_counter == 0) {
LogInfoWithoutPrefix(absl::StrCat("I ", IterationStatsLabelString()));
}
@@ -1081,14 +1085,14 @@ Solver::UpdateIterationStatsAndCheckTermination(
LogInfoWithoutPrefix(absl::StrCat(
"C ", ToString(stats, params_.termination_criteria(),
original_bound_norms_, POINT_TYPE_CURRENT_ITERATE)));
} else if (VLOG_IS_ON(3)) {
} else if (params_.verbosity_level() >= 3) {
if (log_counter == 0) {
LogInfoWithoutPrefix(IterationStatsLabelString());
}
LogInfoWithoutPrefix(ToString(stats, params_.termination_criteria(),
original_bound_norms_,
POINT_TYPE_AVERAGE_ITERATE));
} else {
} else if (params_.verbosity_level() >= 2) {
if (log_counter == 0) {
LogInfoWithoutPrefix(IterationStatsLabelShortString());
}
@@ -1152,12 +1156,14 @@ void Solver::ApplyRestartChoice(const RestartChoice restart_to_apply) {
case RESTART_CHOICE_NO_RESTART:
return;
case RESTART_CHOICE_WEIGHTED_AVERAGE_RESET:
VLOG(4) << "Restarted to current on iteration " << iterations_completed_
<< " after " << primal_average_.NumTerms() << " iterations";
LOG_IF(INFO, params_.verbosity_level() >= 4)
<< "Restarted to current on iteration " << iterations_completed_
<< " after " << primal_average_.NumTerms() << " iterations";
break;
case RESTART_CHOICE_RESTART_TO_AVERAGE:
VLOG(4) << "Restarted to average on iteration " << iterations_completed_
<< " after " << primal_average_.NumTerms() << " iterations";
LOG_IF(INFO, params_.verbosity_level() >= 4)
<< "Restarted to average on iteration " << iterations_completed_
<< " after " << primal_average_.NumTerms() << " iterations";
current_primal_solution_ = primal_average_.ComputeAverage();
current_dual_solution_ = dual_average_.ComputeAverage();
current_dual_product_ = TransposedMatrixVectorProduct(
@@ -1258,52 +1264,53 @@ void Solver::LogInnerIterationLimitHit() const {
}
void Solver::LogQuadraticProgramStats(const QuadraticProgramStats& stats) {
VLOG(1) << absl::StrFormat("There are %i variables, %i constraints, and %i ",
stats.num_variables(), stats.num_constraints(),
stats.constraint_matrix_num_nonzeros())
<< "constraint matrix nonzeros.";
LOG(INFO) << absl::StrFormat(
"There are %i variables, %i constraints, and %i ",
stats.num_variables(), stats.num_constraints(),
stats.constraint_matrix_num_nonzeros())
<< "constraint matrix nonzeros.";
if (WorkingQp().constraint_matrix.nonZeros() > 0) {
VLOG(1) << "Absolute values of nonzero constraint matrix elements: "
<< absl::StrFormat("largest=%f, smallest=%f, avg=%f",
stats.constraint_matrix_abs_max(),
stats.constraint_matrix_abs_min(),
stats.constraint_matrix_abs_avg());
VLOG(1) << "Constraint matrix, infinity norm: "
<< absl::StrFormat("max(row & col)=%f, min_col=%f, min_row=%f",
stats.constraint_matrix_abs_max(),
stats.constraint_matrix_col_min_l_inf_norm(),
stats.constraint_matrix_row_min_l_inf_norm());
VLOG(1) << "Constraint bounds statistics (max absolute value per row): "
<< absl::StrFormat("largest=%f, smallest=%f, avg=%f, l2_norm=%f",
stats.combined_bounds_max(),
stats.combined_bounds_min(),
stats.combined_bounds_avg(),
stats.combined_bounds_l2_norm());
LOG(INFO) << "Absolute values of nonzero constraint matrix elements: "
<< absl::StrFormat("largest=%f, smallest=%f, avg=%f",
stats.constraint_matrix_abs_max(),
stats.constraint_matrix_abs_min(),
stats.constraint_matrix_abs_avg());
LOG(INFO) << "Constraint matrix, infinity norm: "
<< absl::StrFormat("max(row & col)=%f, min_col=%f, min_row=%f",
stats.constraint_matrix_abs_max(),
stats.constraint_matrix_col_min_l_inf_norm(),
stats.constraint_matrix_row_min_l_inf_norm());
LOG(INFO) << "Constraint bounds statistics (max absolute value per row): "
<< absl::StrFormat("largest=%f, smallest=%f, avg=%f, l2_norm=%f",
stats.combined_bounds_max(),
stats.combined_bounds_min(),
stats.combined_bounds_avg(),
stats.combined_bounds_l2_norm());
}
if (!IsLinearProgram(WorkingQp())) {
VLOG(1) << absl::StrFormat(
LOG(INFO) << absl::StrFormat(
"There are %i nonzero diagonal coefficients in the objective matrix.",
stats.objective_matrix_num_nonzeros());
VLOG(1) << "Absolute values of nonzero objective matrix elements: "
<< absl::StrFormat("largest=%f, smallest=%f, avg=%f",
stats.objective_matrix_abs_max(),
stats.objective_matrix_abs_min(),
stats.objective_matrix_abs_avg());
LOG(INFO) << "Absolute values of nonzero objective matrix elements: "
<< absl::StrFormat("largest=%f, smallest=%f, avg=%f",
stats.objective_matrix_abs_max(),
stats.objective_matrix_abs_min(),
stats.objective_matrix_abs_avg());
}
VLOG(1) << "Absolute values of objective vector elements: "
<< absl::StrFormat("largest=%f, smallest=%f, avg=%f, l2_norm=%f",
stats.objective_vector_abs_max(),
stats.objective_vector_abs_min(),
stats.objective_vector_abs_avg(),
stats.objective_vector_l2_norm());
LOG(INFO) << "Absolute values of objective vector elements: "
<< absl::StrFormat("largest=%f, smallest=%f, avg=%f, l2_norm=%f",
stats.objective_vector_abs_max(),
stats.objective_vector_abs_min(),
stats.objective_vector_abs_avg(),
stats.objective_vector_l2_norm());
VLOG(1) << "Gaps between variable upper and lower bounds: "
<< absl::StrFormat(
"#finite=%i of %i, largest=%f, smallest=%f, avg=%f",
stats.variable_bound_gaps_num_finite(), stats.num_variables(),
stats.variable_bound_gaps_max(),
stats.variable_bound_gaps_min(),
stats.variable_bound_gaps_avg());
LOG(INFO) << "Gaps between variable upper and lower bounds: "
<< absl::StrFormat(
"#finite=%i of %i, largest=%f, smallest=%f, avg=%f",
stats.variable_bound_gaps_num_finite(),
stats.num_variables(), stats.variable_bound_gaps_max(),
stats.variable_bound_gaps_min(),
stats.variable_bound_gaps_avg());
}
InnerStepOutcome Solver::TakeMalitskyPockStep() {
@@ -1667,6 +1674,7 @@ PrimalAndDualSolution Solver::RecoverOriginalSolution(
SolverResult Solver::Solve(
absl::optional<PrimalAndDualSolution> initial_solution,
const std::atomic<bool>* interrupt_solve,
IterationStatsCallback iteration_stats_callback) {
SolveLog solve_log;
if (sharded_working_qp_.Qp().problem_name.has_value()) {
@@ -1680,8 +1688,10 @@ SolverResult Solver::Solve(
const std::string preprocessing_string = absl::StrCat(
params_.presolve_options().use_glop() ? "presolving and " : "",
"rescaling:");
VLOG(1) << "Problem stats before " << preprocessing_string;
LogQuadraticProgramStats(solve_log.original_problem_stats());
if (params_.verbosity_level() >= 1) {
LOG(INFO) << "Problem stats before " << preprocessing_string;
LogQuadraticProgramStats(solve_log.original_problem_stats());
}
timer_.Start();
iteration_stats_callback_ = std::move(iteration_stats_callback);
absl::optional<TerminationReason> maybe_terminate =
@@ -1748,8 +1758,10 @@ SolverResult Solver::Solve(
ComputeAndApplyRescaling();
*solve_log.mutable_preprocessed_problem_stats() = ComputeStats(
sharded_working_qp_, params_.infinite_constraint_bound_threshold());
VLOG(1) << "Problem stats after " << preprocessing_string;
LogQuadraticProgramStats(solve_log.preprocessed_problem_stats());
if (params_.verbosity_level() >= 1) {
LOG(INFO) << "Problem stats after " << preprocessing_string;
LogQuadraticProgramStats(solve_log.preprocessed_problem_stats());
}
if (params_.linesearch_rule() ==
PrimalDualHybridGradientParams::CONSTANT_STEP_SIZE_RULE) {
@@ -1817,6 +1829,13 @@ SolverResult Solver::Solve(
if (maybe_result.has_value()) {
return maybe_result.value();
}
// Check for interrupt on every iteration.
if (interrupt_solve != nullptr && interrupt_solve->load() == true) {
return ConstructSolverResult(
PrimalAverage(), DualAverage(),
CreateSimpleIterationStats(RESTART_CHOICE_NO_RESTART),
TERMINATION_REASON_INTERRUPTED_BY_USER, POINT_TYPE_NONE, solve_log);
}
// TODO(user): If we use a step rule that could reject many steps in a
// row, we should add a termination check within this loop also. For the
@@ -1847,14 +1866,17 @@ SolverResult Solver::Solve(
SolverResult PrimalDualHybridGradient(
QuadraticProgram qp, const PrimalDualHybridGradientParams& params,
const std::atomic<bool>* interrupt_solve,
IterationStatsCallback iteration_stats_callback) {
return PrimalDualHybridGradient(std::move(qp), params, absl::nullopt,
interrupt_solve,
std::move(iteration_stats_callback));
}
SolverResult PrimalDualHybridGradient(
QuadraticProgram qp, const PrimalDualHybridGradientParams& params,
absl::optional<PrimalAndDualSolution> initial_solution,
const std::atomic<bool>* interrupt_solve,
IterationStatsCallback iteration_stats_callback) {
const absl::Status params_status =
ValidatePrimalDualHybridGradientParams(params);
@@ -1881,7 +1903,7 @@ SolverResult PrimalDualHybridGradient(
"The objective scaling factor cannot be zero.");
}
Solver solver(std::move(qp), params);
return solver.Solve(std::move(initial_solution),
return solver.Solve(std::move(initial_solution), interrupt_solve,
std::move(iteration_stats_callback));
}

17
ortools/pdlp/primal_dual_hybrid_gradient.h
View File

@@ -14,6 +14,7 @@
#ifndef PDLP_PRIMAL_DUAL_HYBRID_GRADIENT_H_
#define PDLP_PRIMAL_DUAL_HYBRID_GRADIENT_H_
#include <atomic>
#include <functional>
#include "Eigen/Core"
@@ -50,10 +51,10 @@ struct PrimalAndDualSolution {
// infeasible but no certificate of infeasibility is available. The
// primal_solution and dual_solution have no meaning. This status is only used
// when presolve is enabled.
// * TIME_LIMIT, ITERATION_LIMIT, KKT_MATRIX_PASS_LIMIT, NUMERICAL_ERROR:
// the vectors contain an iterate at the time that the respective event
// occurred. Their values may or may not be meaningful. In some cases
// solution quality information is available; see documentation for
// * TIME_LIMIT, ITERATION_LIMIT, KKT_MATRIX_PASS_LIMIT, NUMERICAL_ERROR,
// INTERRUPTED_BY_USER: the vectors contain an iterate at the time that the
// respective event occurred. Their values may or may not be meaningful. In
// some cases solution quality information is available; see documentation for
// solve_log.solution_type.
// * INVALID_PROBLEM, INVALID_PARAMETER, OTHER: the solution vectors are
// meaningless and may not have lengths consistent with the input problem.
@@ -100,8 +101,12 @@ struct IterationCallbackInfo {
// and params.primal_weight_update_smoothing controls how primal_weight is
// updated.
//
// If interrupt_solve is not nullptr, then the solver will periodically check if
// interrupt_solve->load() is true, in which case the solve will terminate with
// TERMINATION_REASON_INTERRUPTED_BY_USER.
//
// If iteration_stats_callback is not nullptr, then at each termination step
// (when iteration stats are VLOG(2)ed), iteration_stats_callback will also
// (when iteration stats are logged), iteration_stats_callback will also
// be called with those iteration stats.
//
// Callers MUST check solve_log.termination_reason before using the vectors in
@@ -115,6 +120,7 @@ struct IterationCallbackInfo {
// modifies its copy.
SolverResult PrimalDualHybridGradient(
QuadraticProgram qp, const PrimalDualHybridGradientParams& params,
const std::atomic<bool>* interrupt_solve = nullptr,
std::function<void(const IterationCallbackInfo&)> iteration_stats_callback =
nullptr);
@@ -126,6 +132,7 @@ SolverResult PrimalDualHybridGradient(
SolverResult PrimalDualHybridGradient(
QuadraticProgram qp, const PrimalDualHybridGradientParams& params,
absl::optional<PrimalAndDualSolution> initial_solution,
const std::atomic<bool>* interrupt_solve = nullptr,
std::function<void(const IterationCallbackInfo&)> iteration_stats_callback =
nullptr);

50
ortools/pdlp/primal_dual_hybrid_gradient_test.cc
View File

@@ -14,6 +14,7 @@
#include "ortools/pdlp/primal_dual_hybrid_gradient.h"
#include <algorithm>
#include <atomic>
#include <cmath>
#include <cstdint>
#include <limits>
@@ -1016,7 +1017,7 @@ TEST(PrimalDualHybridGradientTest, CallsCallback) {
int callback_count = 0;
SolverResult output = PrimalDualHybridGradient(
TestLp(), params,
TestLp(), params, /*interrupt_solve=*/nullptr,
[&callback_count](const IterationCallbackInfo& callback_info) {
++callback_count;
});
@@ -1082,6 +1083,53 @@ TEST(PrimalDualHybridGradientTest, EmptyQp) {
EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}
TEST(PrimalDualHybridGradientTest, RespectsInterrupt) {
std::atomic<bool> interrupt_solve;
PrimalDualHybridGradientParams params;
params.mutable_termination_criteria()->set_eps_optimal_absolute(0.0);
params.mutable_termination_criteria()->set_eps_optimal_relative(0.0);
interrupt_solve.store(true);
const SolverResult output =
PrimalDualHybridGradient(TestLp(), params, &interrupt_solve);
EXPECT_EQ(output.solve_log.termination_reason(),
TERMINATION_REASON_INTERRUPTED_BY_USER);
}
TEST(PrimalDualHybridGradientTest, RespectsInterruptFromCallback) {
std::atomic<bool> interrupt_solve;
PrimalDualHybridGradientParams params;
params.mutable_termination_criteria()->set_eps_optimal_absolute(0.0);
params.mutable_termination_criteria()->set_eps_optimal_relative(0.0);
interrupt_solve.store(false);
auto callback = [&](const IterationCallbackInfo& info) {
if (info.iteration_stats.iteration_number() >= 10) {
interrupt_solve.store(true);
}
};
const SolverResult output =
PrimalDualHybridGradient(TestLp(), params, &interrupt_solve, callback);
EXPECT_EQ(output.solve_log.termination_reason(),
TERMINATION_REASON_INTERRUPTED_BY_USER);
EXPECT_GE(output.solve_log.iteration_count(), 10);
}
TEST(PrimalDualHybridGradientTest, IgnoresFalseInterrupt) {
std::atomic<bool> interrupt_solve;
PrimalDualHybridGradientParams params;
params.mutable_termination_criteria()->set_eps_optimal_absolute(0.0);
params.mutable_termination_criteria()->set_eps_optimal_relative(0.0);
params.mutable_termination_criteria()->set_kkt_matrix_pass_limit(1);
interrupt_solve.store(false);
const SolverResult output =
PrimalDualHybridGradient(TestLp(), params, &interrupt_solve);
EXPECT_EQ(output.solve_log.termination_reason(),
TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT);
}
// Verifies that the primal and dual solution satisfy the bounds constraints.
// This function uses ASSERTS rather than EXPECTs because it's used with large
// QPs that would spam the logs otherwise.

1
ortools/pdlp/solve_log.proto
View File

@@ -306,6 +306,7 @@ enum TerminationReason {
TERMINATION_REASON_TIME_LIMIT = 4;
TERMINATION_REASON_ITERATION_LIMIT = 5;
TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT = 8;
TERMINATION_REASON_INTERRUPTED_BY_USER = 12;
TERMINATION_REASON_NUMERICAL_ERROR = 6;
// Indicates that the solver detected invalid problem data, e.g., inconsistent
// bounds.

13
ortools/pdlp/solvers.proto
View File

@@ -210,6 +210,19 @@ message PrimalDualHybridGradientParams {
// increase the size of the output.
optional bool record_iteration_stats = 3;
// The verbosity of logging.
// 0: No informational logging. (Errors are logged.)
// 1: Summary statistics only. No iteration-level details.
// 2: A table of iteration-level statistics is logged.
// (See ToShortString() in primal_dual_hybrid_gradient.cc).
// 3: A more detailed table of iteration-level statistics is logged.
// (See ToString() in primal_dual_hybrid_gradient.cc).
// 4: For iteration-level details, prints the statistics of both the average
// (prefixed with A) and the current iterate (prefixed with C). Also prints
// internal algorithmic state and details.
// Logging at levels 2-4 also includes messages from level 1.
optional int32 verbosity_level = 26 [default = 0];
// The frequency at which extra work is performed to make major algorithmic
// decisions, e.g., performing restarts and updating the primal weight. Major
// iterations also trigger a termination check. For best performance using

3
ortools/pdlp/solvers_proto_validation.cc
View File

@@ -87,6 +87,9 @@ absl::Status ValidatePrimalDualHybridGradientParams(
if (params.num_threads() <= 0) {
return InvalidArgumentError("num_threads must be positive");
}
if (params.verbosity_level() < 0) {
return InvalidArgumentError("verbosity_level must be non-negative");
}
if (params.major_iteration_frequency() <= 0) {
return InvalidArgumentError("major_iteration_frequency must be positive");
}

8
ortools/pdlp/solvers_proto_validation_test.cc
View File

@@ -179,6 +179,14 @@ TEST(ValidatePrimalDualHybridGradientParams, BadNumThreads) {
EXPECT_THAT(status.message(), HasSubstr("num_threads"));
}
TEST(ValidatePrimalDualHybridGradientParams, BadVerbosityLevel) {
PrimalDualHybridGradientParams params;
params.set_verbosity_level(-1);
const absl::Status status = ValidatePrimalDualHybridGradientParams(params);
EXPECT_EQ(status.code(), absl::StatusCode::kInvalidArgument);
EXPECT_THAT(status.message(), HasSubstr("verbosity_level"));
}
TEST(ValidatePrimalDualHybridGradientParams, BadMajorIterationFrequency) {
PrimalDualHybridGradientParams params;
params.set_major_iteration_frequency(0);

22
ortools/pdlp/trust_region_test.cc
View File

@@ -651,7 +651,7 @@ TEST_P(ComputeLocalizedLagrangianBoundsTest, OptimalInBoundRange) {
const double primal_distance_squared_to_optimal =
0.5 * (1.0 + 8.0 * 8.0 + 1.0 + 0.5 * 0.5);
const double dual_distance_squared_to_optimal =
0.5 * (4.0 + 2.375 * 2.375 + 4.0 / 6.0);
0.5 * (4.0 + 2.375 * 2.375 + 4.0 / 9.0);
const double distance_to_optimal =
primal_dual_norm == PrimalDualNorm::kEuclideanNorm
? std::sqrt(primal_distance_squared_to_optimal +
@@ -774,19 +774,23 @@ TEST_P(ComputeLocalizedLagrangianBoundsTest, AcceptsCachedProducts) {
primal_product << 6.0, 0.0, 0.0, -3.0;
VectorXd dual_product = VectorXd::Zero(4);
const double primal_distance_to_optimal =
std::sqrt(0.5 * (1.0 + 8.0 * 8.0 + 1.0 + 0.5 * 0.5));
const double dual_distance_to_optimal =
std::sqrt(0.5 * (4.0 + 2.375 * 2.375 + 4.0 / 6.0));
const double max_norm_distance_to_optimal =
std::max(primal_distance_to_optimal, dual_distance_to_optimal);
const auto [primal_dual_norm, use_diagonal_qp_solver] = GetParam();
const double primal_distance_squared_to_optimal =
0.5 * (1.0 + 8.0 * 8.0 + 1.0 + 0.5 * 0.5);
const double dual_distance_squared_to_optimal =
0.5 * (4.0 + 2.375 * 2.375 + 4.0 / 9.0);
const double distance_to_optimal =
primal_dual_norm == PrimalDualNorm::kEuclideanNorm
? std::sqrt(primal_distance_squared_to_optimal +
dual_distance_squared_to_optimal)
: std::sqrt(std::max(primal_distance_squared_to_optimal,
dual_distance_squared_to_optimal));
LocalizedLagrangianBounds bounds = ComputeLocalizedLagrangianBounds(
lp, primal_solution, dual_solution, primal_dual_norm,
/*primal_weight=*/1.0,
/*radius=*/max_norm_distance_to_optimal,
/*radius=*/distance_to_optimal,
/*primal_product=*/&primal_product,
/*dual_product=*/&dual_product, use_diagonal_qp_solver,
/*diagonal_qp_trust_region_solver_tolerance=*/1.0e-6);