diff --git a/makefiles/Makefile.dotnet.mk b/makefiles/Makefile.dotnet.mk index 60b5f01173..565e2c4808 100644 --- a/makefiles/Makefile.dotnet.mk +++ b/makefiles/Makefile.dotnet.mk @@ -703,6 +703,8 @@ clean_dotnet: -$(DELREC) ortools$Slinear_solver$Ssamples$Sobj -$(DELREC) ortools$Ssat$Ssamples$Sbin -$(DELREC) ortools$Ssat$Ssamples$Sobj + -$(DEL) *.nupkg + -$(DEL) *.snupkg -@"$(DOTNET_BIN)" nuget locals all --clear ############# diff --git a/makefiles/Makefile.java.mk b/makefiles/Makefile.java.mk index 9e3ea45a0c..0b16550296 100644 --- a/makefiles/Makefile.java.mk +++ b/makefiles/Makefile.java.mk @@ -494,6 +494,7 @@ clean_java: -$(DEL) $(OBJ_DIR)$Sswig$S*_java_wrap.$O -$(DEL) $(LIB_DIR)$S$(LIB_PREFIX)jni*.$(JNI_LIB_EXT) -$(DEL) $(LIB_DIR)$S*.jar + -$(DEL) *.jar -$(DELREC) $(TEMP_JAVA_DIR) ############# diff --git a/makefiles/Makefile.python.mk b/makefiles/Makefile.python.mk index f7fe5ede9c..3b06f2eb7f 100644 --- a/makefiles/Makefile.python.mk +++ b/makefiles/Makefile.python.mk @@ -501,7 +501,7 @@ test_package_python: package_python $(COPY) ortools$Sconstraint_solver$Ssamples$Scvrptw_break.py $(PYPI_ARCHIVE_TEMP_DIR)$Svenv ifneq ($(SYSTEM),win) $(PYPI_ARCHIVE_TEMP_DIR)/venv/bin/python -m pip install $(PYPI_ARCHIVE_TEMP_DIR)/ortools/dist/*.whl - $(PYPI_ARCHIVE_TEMP_DIR)/venv/bin/python -m pip install pandas matplotlib + $(PYPI_ARCHIVE_TEMP_DIR)/venv/bin/python -m pip install pandas matplotlibgit $(PYPI_ARCHIVE_TEMP_DIR)/venv/bin/python $(PYPI_ARCHIVE_TEMP_DIR)/venv/test.py $(PYPI_ARCHIVE_TEMP_DIR)/venv/bin/python $(PYPI_ARCHIVE_TEMP_DIR)/venv/simple_knapsack_program.py $(PYPI_ARCHIVE_TEMP_DIR)/venv/bin/python $(PYPI_ARCHIVE_TEMP_DIR)/venv/simple_max_flow_program.py @@ -678,6 +678,7 @@ clean_python: -$(DEL) $(GEN_PATH)$Sortools$Sutil$S_* -$(DEL) $(LIB_DIR)$S_*.$(SWIG_PYTHON_LIB_SUFFIX) -$(DEL) $(OBJ_DIR)$Sswig$S*python_wrap.$O + -$(DEL) *.whl -$(DELREC) temp_python* ############# diff --git a/ortools/linear_solver/linear_solver.proto b/ortools/linear_solver/linear_solver.proto index 1aacecec98..d6d1716ed4 100644 --- a/ortools/linear_solver/linear_solver.proto +++ b/ortools/linear_solver/linear_solver.proto @@ -603,7 +603,7 @@ message MPSolveInfo { optional double solve_user_time_seconds = 2; } -// Next id: 11. +// Next id: 12. message MPSolutionResponse { // Result of the optimization. optional /*required*/ MPSolverResponseStatus status = 1 @@ -637,6 +637,10 @@ message MPSolutionResponse { // by SCIP and Gurobi proto solves. optional MPSolveInfo solve_info = 10; + // Opaque solver-specific information. + // For the PDLP solver, this is a serialized pdlp::SolveLog proto. + optional string solver_specific_info = 11; + // [Advanced usage.] // Values of the dual variables values in the same order as the // MPModelProto::constraint field. This is a dense representation. diff --git a/ortools/linear_solver/pdlp_interface.cc b/ortools/linear_solver/pdlp_interface.cc index 16f8bf4022..99e5b007c2 100644 --- a/ortools/linear_solver/pdlp_interface.cc +++ b/ortools/linear_solver/pdlp_interface.cc @@ -129,25 +129,27 @@ MPSolver::ResultStatus PdlpInterface::Solve(const MPSolverParameters& param) { MPModelProto model_proto; solver_->ExportModelToProto(&model_proto); // TODO(user): Adjust logging level based on quiet_. - absl::StatusOr solution = + absl::StatusOr response = pdlp::SolveMpModelProto(model_proto, parameters_, /*relax_integer_variables=*/true); - if (!solution.ok()) { - LOG(ERROR) << "Unexpected error solving with PDLP: " << solution.status(); + if (!response.ok()) { + LOG(ERROR) << "Unexpected error solving with PDLP: " << response.status(); return MPSolver::ABNORMAL; } - const MPSolutionResponse& response = solution->response; - // The solution must be marked as synchronized even when no solution exists. sync_status_ = SOLUTION_SYNCHRONIZED; - result_status_ = static_cast(response.status()); - solve_log_ = solution->solve_log; + result_status_ = static_cast(response->status()); + LOG_IF(DFATAL, !response->has_solver_specific_info()) + << response->DebugString(); + if (!solve_log_.ParseFromString(response->solver_specific_info())) { + LOG(DFATAL) << "Unable to parse PDLP's SolveLog from solver_specific_info"; + } - if (response.status() == MPSOLVER_FEASIBLE || - response.status() == MPSOLVER_OPTIMAL) { - const absl::Status result = solver_->LoadSolutionFromProto(response); + if (response->status() == MPSOLVER_FEASIBLE || + response->status() == MPSOLVER_OPTIMAL) { + const absl::Status result = solver_->LoadSolutionFromProto(*response); if (!result.ok()) { LOG(ERROR) << "LoadSolutionFromProto failed: " << result; } @@ -162,13 +164,13 @@ absl::optional PdlpInterface::DirectlySolveProto( return absl::nullopt; } - MPSolutionResponse response; + MPSolutionResponse error_response; if (request.has_solver_specific_parameters()) { if (!solver_->SetSolverSpecificParametersAsString( request.solver_specific_parameters())) { - response.set_status( + error_response.set_status( MPSolverResponseStatus::MPSOLVER_MODEL_INVALID_SOLVER_PARAMETERS); - return response; + return error_response; } } @@ -182,27 +184,28 @@ absl::optional PdlpInterface::DirectlySolveProto( } const absl::optional> optional_model = - ExtractValidMPModelOrPopulateResponseStatus(request, &response); + ExtractValidMPModelOrPopulateResponseStatus(request, &error_response); if (!optional_model) { LOG_IF(WARNING, request.enable_internal_solver_output()) << "Failed to extract a valid model from protocol buffer. Status: " - << ProtoEnumToString(response.status()) << " (" - << response.status() << "): " << response.status_str(); + << ProtoEnumToString(error_response.status()) + << " (" << error_response.status() + << "): " << error_response.status_str(); return absl::nullopt; } - absl::StatusOr solution = + absl::StatusOr response = pdlp::SolveMpModelProto(optional_model->get(), parameters_, /*relax_integer_variables=*/true); - if (!solution.ok()) { - LOG(ERROR) << "Unexpected error solving with PDLP: " << solution.status(); - response.set_status(MPSolverResponseStatus::MPSOLVER_ABNORMAL); - response.set_status_str(solution.status().ToString()); - return response; + if (!response.ok()) { + LOG(ERROR) << "Unexpected error solving with PDLP: " << response.status(); + error_response.set_status(MPSolverResponseStatus::MPSOLVER_ABNORMAL); + error_response.set_status_str(response.status().ToString()); + return error_response; } - return solution->response; + return *response; } void PdlpInterface::Reset() { ResetExtractionInformation(); } diff --git a/ortools/pdlp/BUILD.bazel b/ortools/pdlp/BUILD.bazel index 43b4ff165f..327be3893b 100644 --- a/ortools/pdlp/BUILD.bazel +++ b/ortools/pdlp/BUILD.bazel @@ -1,16 +1,3 @@ -# Copyright 2010-2021 Google LLC -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. - load("@rules_cc//cc:defs.bzl", "cc_proto_library") package(default_visibility = ["//visibility:public"]) @@ -125,7 +112,6 @@ cc_library( cc_test( name = "test_util_test", srcs = ["test_util_test.cc"], - defines = ["_USE_MATH_DEFINES"], deps = [ ":gtest_main", ":test_util", @@ -133,6 +119,7 @@ cc_test( "@com_google_absl//absl/types:span", "@eigen//:eigen3", ], + defines = ["_USE_MATH_DEFINES"], ) cc_library( @@ -301,6 +288,7 @@ cc_library( ], ) + cc_library( name = "iteration_stats", srcs = ["iteration_stats.cc"], diff --git a/ortools/pdlp/primal_dual_hybrid_gradient.cc b/ortools/pdlp/primal_dual_hybrid_gradient.cc index b36d6bab3f..76f1ea3a88 100644 --- a/ortools/pdlp/primal_dual_hybrid_gradient.cc +++ b/ortools/pdlp/primal_dual_hybrid_gradient.cc @@ -1886,7 +1886,7 @@ SolverResult PrimalDualHybridGradient( std::move(iteration_stats_callback)); } -absl::StatusOr SolveMpModelProto( +absl::StatusOr SolveMpModelProto( const MPModelProto& model, const PrimalDualHybridGradientParams& params, const bool relax_integer_variables, IterationStatsCallback iteration_stats_callback) { @@ -1944,8 +1944,9 @@ absl::StatusOr SolveMpModelProto( response.add_reduced_cost(objective_scaling_factor * v); } - return SolutionResponseAndLog{.response = response, - .solve_log = pdhg_result.solve_log}; + response.set_solver_specific_info(pdhg_result.solve_log.SerializeAsString()); + + return response; } namespace internal { diff --git a/ortools/pdlp/primal_dual_hybrid_gradient.h b/ortools/pdlp/primal_dual_hybrid_gradient.h index 27d301ec5c..e5b0ab3602 100644 --- a/ortools/pdlp/primal_dual_hybrid_gradient.h +++ b/ortools/pdlp/primal_dual_hybrid_gradient.h @@ -129,22 +129,16 @@ SolverResult PrimalDualHybridGradient( std::function iteration_stats_callback = nullptr); -struct SolutionResponseAndLog { - MPSolutionResponse response; - SolveLog solve_log; -}; - // Uses PrimalDualHybridGradient to solve the instance specified by the // MPModelProto. If relax_integer_variables is true, integrality constraints are // relaxed before solving. If false, integrality constraints result in an error. -// The MPSolutionResponse in the result enables partial compatibility with the -// MPSolver ecosystem, e.g., with MPSolver::LoadSolutionFromProto. The SolveLog -// has more detailed solve information. 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 QuadraticProgram fails. The lack of an error does not imply -// success. Check solve_log.termination_reason to see if the solve succeeded. -absl::StatusOr SolveMpModelProto( +// 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 QuadraticProgram fails. The +// lack of an error does not imply success. Check the SolveLog's +// termination_reason for more refined status details. +absl::StatusOr SolveMpModelProto( const MPModelProto& model, const PrimalDualHybridGradientParams& params, bool relax_integer_variables = false, std::function iteration_stats_callback = diff --git a/ortools/pdlp/primal_dual_hybrid_gradient_test.cc b/ortools/pdlp/primal_dual_hybrid_gradient_test.cc index 690f0cc839..af17e49768 100644 --- a/ortools/pdlp/primal_dual_hybrid_gradient_test.cc +++ b/ortools/pdlp/primal_dual_hybrid_gradient_test.cc @@ -1126,24 +1126,26 @@ TEST(SolveProtoTest, SolvesEasyMinimizationLp) { )pb"); // Using default convergence tolerances and 1 thread. PrimalDualHybridGradientParams params; - absl::StatusOr result = + absl::StatusOr response = SolveMpModelProto(proto, params); - ASSERT_TRUE(result.ok()) << result.status(); - EXPECT_EQ(result->response.status(), MPSOLVER_OPTIMAL); - EXPECT_DOUBLE_EQ(result->response.objective_value(), -2.0); - EXPECT_THAT(result->response.variable_value(), ElementsAre(0.0, 1.0)); - EXPECT_THAT(result->response.dual_value(), ElementsAre(-2.0)); - EXPECT_THAT(result->response.reduced_cost(), ElementsAre(2.0, 0.0)); + ASSERT_TRUE(response.ok()) << response.status(); + SolveLog solve_log; + ASSERT_TRUE(solve_log.ParseFromString(response->solver_specific_info())); + EXPECT_EQ(response->status(), MPSOLVER_OPTIMAL); + EXPECT_DOUBLE_EQ(response->objective_value(), -2.0); + EXPECT_THAT(response->variable_value(), ElementsAre(0.0, 1.0)); + EXPECT_THAT(response->dual_value(), ElementsAre(-2.0)); + EXPECT_THAT(response->reduced_cost(), ElementsAre(2.0, 0.0)); - EXPECT_TRUE(result->solve_log.has_solution_stats()); - EXPECT_GT(result->solve_log.iteration_count(), 0); + EXPECT_TRUE(solve_log.has_solution_stats()); + EXPECT_GT(solve_log.iteration_count(), 0); // The signs of the duals are consistent with Glop. const MPSolutionResponse glop_response = GlopSolution(proto); EXPECT_THAT(glop_response.dual_value(), - ElementsAreArray(result->response.dual_value())); + ElementsAreArray(response->dual_value())); EXPECT_THAT(glop_response.reduced_cost(), - ElementsAreArray(result->response.reduced_cost())); + ElementsAreArray(response->reduced_cost())); } // This test solves an LP that's designed to terminate with @@ -1167,16 +1169,19 @@ TEST(SolveProtoTest, SolvesWithDifferenceOfIterates) { params.set_major_iteration_frequency(1); params.set_l_inf_ruiz_iterations(0); params.set_l2_norm_rescaling(false); - absl::StatusOr result = + absl::StatusOr response = SolveMpModelProto(proto, params); - ASSERT_TRUE(result.ok()) << result.status(); - EXPECT_EQ(result->response.status(), MPSOLVER_NOT_SOLVED); - EXPECT_THAT(result->response.variable_value(), ElementsAre(-0.4, 0.0)); - EXPECT_THAT(result->response.dual_value(), ElementsAre(0.0)); - EXPECT_THAT(result->response.reduced_cost(), ElementsAre(0.0, 0.0)); - EXPECT_FALSE(result->response.has_objective_value()); - EXPECT_TRUE(result->solve_log.has_solution_stats()); - EXPECT_EQ(result->solve_log.iteration_count(), 1); + ASSERT_TRUE(response.ok()) << response.status(); + SolveLog solve_log; + ASSERT_TRUE(solve_log.ParseFromString(response->solver_specific_info())); + + EXPECT_EQ(response->status(), MPSOLVER_NOT_SOLVED); + EXPECT_THAT(response->variable_value(), ElementsAre(-0.4, 0.0)); + EXPECT_THAT(response->dual_value(), ElementsAre(0.0)); + EXPECT_THAT(response->reduced_cost(), ElementsAre(0.0, 0.0)); + EXPECT_FALSE(response->has_objective_value()); + EXPECT_TRUE(solve_log.has_solution_stats()); + EXPECT_EQ(solve_log.iteration_count(), 1); } TEST_P(SolveMpModelProtoMaybePresolveTest, SolvesEasyMaximizationLp) { @@ -1201,26 +1206,30 @@ TEST_P(SolveMpModelProtoMaybePresolveTest, SolvesEasyMaximizationLp) { const bool presolve_on = GetParam(); params.mutable_presolve_options()->set_use_glop(presolve_on); - absl::StatusOr result = + absl::StatusOr response = SolveMpModelProto(proto, params); - ASSERT_TRUE(result.ok()) << result.status(); - EXPECT_EQ(result->response.status(), MPSOLVER_OPTIMAL); - EXPECT_DOUBLE_EQ(result->response.objective_value(), 2.0); - EXPECT_THAT(result->response.variable_value(), ElementsAre(0.0, 1.0)); - EXPECT_THAT(result->response.dual_value(), ElementsAre(2.0)); - EXPECT_THAT(result->response.reduced_cost(), ElementsAre(-2.0, 0.0)); + ASSERT_TRUE(response.ok()) << response.status(); + SolveLog solve_log; + ASSERT_TRUE(solve_log.ParseFromString(response->solver_specific_info())); - EXPECT_TRUE(result->solve_log.has_solution_stats()); + EXPECT_EQ(response->status(), MPSOLVER_OPTIMAL); + EXPECT_DOUBLE_EQ(response->objective_value(), 2.0); + EXPECT_THAT(response->variable_value(), ElementsAre(0.0, 1.0)); + EXPECT_THAT(response->dual_value(), ElementsAre(2.0)); + EXPECT_THAT(response->reduced_cost(), ElementsAre(-2.0, 0.0)); + + EXPECT_TRUE(solve_log.has_solution_stats()); if (presolve_on) { - EXPECT_EQ(result->solve_log.iteration_count(), 0); + EXPECT_EQ(solve_log.iteration_count(), 0); } else { - EXPECT_GT(result->solve_log.iteration_count(), 0); + EXPECT_GT(solve_log.iteration_count(), 0); } // The signs of the duals are consistent with Glop. - EXPECT_THAT(result->response.dual_value(), - ElementsAreArray(result->response.dual_value())); - EXPECT_THAT(result->response.reduced_cost(), - ElementsAreArray(result->response.reduced_cost())); + const MPSolutionResponse glop_response = GlopSolution(proto); + EXPECT_THAT(glop_response.dual_value(), + ElementsAreArray(response->dual_value())); + EXPECT_THAT(glop_response.reduced_cost(), + ElementsAreArray(response->reduced_cost())); } TEST(SolveProtoTest, RelaxesIntegerVariables) { @@ -1237,13 +1246,13 @@ TEST(SolveProtoTest, RelaxesIntegerVariables) { )pb"); // Using default convergence tolerances and 1 thread. PrimalDualHybridGradientParams params; - absl::StatusOr result = + absl::StatusOr response = SolveMpModelProto(proto, params, /*relax_integer_variables=*/true); - ASSERT_TRUE(result.ok()) << result.status(); - EXPECT_EQ(result->response.status(), MPSOLVER_OPTIMAL); - EXPECT_DOUBLE_EQ(result->response.objective_value(), 1.0); - EXPECT_THAT(result->response.variable_value(), ElementsAre(1.0)); - EXPECT_THAT(result->response.reduced_cost(), ElementsAre(1.0)); + ASSERT_TRUE(response.ok()) << response.status(); + EXPECT_EQ(response->status(), MPSOLVER_OPTIMAL); + EXPECT_DOUBLE_EQ(response->objective_value(), 1.0); + EXPECT_THAT(response->variable_value(), ElementsAre(1.0)); + EXPECT_THAT(response->reduced_cost(), ElementsAre(1.0)); } TEST(SolveProtoTest, ErrorsOnIntegerVariables) { @@ -1289,11 +1298,10 @@ TEST(SolveProtoTest, SolvesEasyLpWithMPSolver) { // Using default convergence tolerances and 1 thread. PrimalDualHybridGradientParams params; - absl::StatusOr result = + absl::StatusOr response = SolveMpModelProto(proto, params); - ASSERT_TRUE(result.ok()) << result.status(); - const absl::Status load_status = - model.LoadSolutionFromProto(result->response); + ASSERT_TRUE(response.ok()) << response.status(); + const absl::Status load_status = model.LoadSolutionFromProto(*response); ASSERT_TRUE(load_status.ok()) << load_status; EXPECT_DOUBLE_EQ(x->solution_value(), 0.0); diff --git a/ortools/pdlp/quadratic_program_io.cc b/ortools/pdlp/quadratic_program_io.cc index e86472af68..fdb697cf5e 100644 --- a/ortools/pdlp/quadratic_program_io.cc +++ b/ortools/pdlp/quadratic_program_io.cc @@ -31,7 +31,6 @@ #include "absl/types/optional.h" #include "ortools/base/basictypes.h" #include "ortools/base/file.h" -#include "ortools/base/helpers.h" #include "ortools/base/logging.h" #include "ortools/base/mathutil.h" #include "ortools/base/status_macros.h" diff --git a/ortools/pdlp/sharded_optimization_utils.cc b/ortools/pdlp/sharded_optimization_utils.cc index f73f91d41b..ccdf15a76d 100644 --- a/ortools/pdlp/sharded_optimization_utils.cc +++ b/ortools/pdlp/sharded_optimization_utils.cc @@ -489,19 +489,6 @@ LagrangianPart ComputePrimalGradient(const ShardedQuadraticProgram& sharded_qp, return result; } -namespace { - -// Returns a subderivative of the concave dual penalty function that appears in -// the Lagrangian: -p(dual; -constraint_upper_bound, -constraint_lower_bound) = -// { constraint_upper_bound * dual when dual < 0, 0 when dual == 0, and -// constraint_lower_bound * dual when dual > 0} -// (as defined at https://developers.google.com/optimization/lp/pdlp_math). -// The subderivative is not necessarily unique when dual == 0. In this case, if -// only one of the bounds is finite, we return that one. If both are finite, we -// return the primal product projected onto the bounds, which causes the dual -// Lagrangian gradient to be zero when the constraint is not violated. If both -// are infinite, we return zero. The value returned is valid only when the -// function is finite-valued. double DualSubgradientCoefficient(const double constraint_lower_bound, const double constraint_upper_bound, const double dual, @@ -528,8 +515,6 @@ double DualSubgradientCoefficient(const double constraint_lower_bound, } } -} // namespace - LagrangianPart ComputeDualGradient(const ShardedQuadraticProgram& sharded_qp, const Eigen::VectorXd& dual_solution, const Eigen::VectorXd& primal_product) { diff --git a/ortools/pdlp/sharded_optimization_utils.h b/ortools/pdlp/sharded_optimization_utils.h index 09384bcbe8..f952d41cd8 100644 --- a/ortools/pdlp/sharded_optimization_utils.h +++ b/ortools/pdlp/sharded_optimization_utils.h @@ -137,6 +137,22 @@ LagrangianPart ComputePrimalGradient(const ShardedQuadraticProgram& sharded_qp, const Eigen::VectorXd& primal_solution, const Eigen::VectorXd& dual_product); +// Returns a subderivative of the concave dual penalty function that appears in +// the Lagrangian: -p(dual; -constraint_upper_bound, -constraint_lower_bound) = +// { constraint_upper_bound * dual when dual < 0, 0 when dual == 0, and +// constraint_lower_bound * dual when dual > 0} +// (as defined at https://developers.google.com/optimization/lp/pdlp_math). +// The subderivative is not necessarily unique when dual == 0. In this case, if +// only one of the bounds is finite, we return that one. If both are finite, we +// return the primal product projected onto the bounds, which causes the dual +// Lagrangian gradient to be zero when the constraint is not violated. If both +// are infinite, we return zero. The value returned is valid only when the +// function is finite-valued. +double DualSubgradientCoefficient(const double constraint_lower_bound, + const double constraint_upper_bound, + const double dual, + const double primal_product); + // Computes the value of the dual part of the Lagrangian function defined at // https://developers.google.com/optimization/lp/pdlp_math, i.e., -h^*(y) and // the gradient of the Lagrangian with respect to the dual variables y, i.e., diff --git a/ortools/pdlp/trust_region.cc b/ortools/pdlp/trust_region.cc index f839c6493c..ed22609332 100644 --- a/ortools/pdlp/trust_region.cc +++ b/ortools/pdlp/trust_region.cc @@ -21,7 +21,6 @@ #include #include "Eigen/Core" -#include "absl/algorithm/container.h" #include "absl/types/optional.h" #include "ortools/base/logging.h" #include "ortools/base/mathutil.h" @@ -98,80 +97,6 @@ class VectorTrustRegionProblem { const VectorXd& norm_weight_; }; -// PrimalTrustRegionProblem defines the primal trust region problem given a -// QuadraticProgram, primal solution, and primal gradient. It captures const -// references to the constructor arguments, which should outlive the class -// instance. -// The corresponding trust region problem is -// min_x primal_gradient' * (x - primal_solution) -// s.t. qp.variable_lower_bounds <= x <= qp.variable_upper_bounds -// || x - primal_solution ||_2 <= target_radius -class PrimalTrustRegionProblem { - public: - PrimalTrustRegionProblem(const QuadraticProgram* qp, - const VectorXd* primal_solution, - const VectorXd* primal_gradient) - : qp_(*qp), - primal_solution_(*primal_solution), - primal_gradient_(*primal_gradient) {} - double Objective(int64_t index) const { return primal_gradient_[index]; } - double LowerBound(int64_t index) const { - return qp_.variable_lower_bounds[index]; - } - double UpperBound(int64_t index) const { - return qp_.variable_upper_bounds[index]; - } - double CenterPoint(int64_t index) const { return primal_solution_[index]; } - double NormWeight(int64_t index) const { return 1.0; } - - private: - const QuadraticProgram& qp_; - const VectorXd& primal_solution_; - const VectorXd& primal_gradient_; -}; - -// DualTrustRegionProblem defines the dual trust region problem given a -// QuadraticProgram, dual solution, and dual gradient. It captures const -// references to the constructor arguments, which should outlive the class -// instance. -// The corresponding trust region problem is -// max_y dual_gradient' * (y - dual_solution) -// s.t. qp.implicit_dual_lower_bounds <= y <= qp.implicit_dual_upper_bounds -// || y - dual_solution ||_2 <= target_radius -// where the implicit dual bounds are those given in -// https://developers.google.com/optimization/lp/pdlp_math#dual_variable_bounds -class DualTrustRegionProblem { - public: - DualTrustRegionProblem(const QuadraticProgram* qp, - const VectorXd* dual_solution, - const VectorXd* dual_gradient) - : qp_(*qp), - dual_solution_(*dual_solution), - dual_gradient_(*dual_gradient) {} - double Objective(int64_t index) const { - // The objective is negated because the trust region problem objective is - // minimize, but for the dual problem we want to maximize the gradient. - return -dual_gradient_[index]; - } - double LowerBound(int64_t index) const { - return std::isfinite(qp_.constraint_upper_bounds[index]) - ? -std::numeric_limits::infinity() - : 0.0; - } - double UpperBound(int64_t index) const { - return std::isfinite(qp_.constraint_lower_bounds[index]) - ? std::numeric_limits::infinity() - : 0.0; - } - double CenterPoint(int64_t index) const { return dual_solution_[index]; } - double NormWeight(int64_t index) const { return 1.0; } - - private: - const QuadraticProgram& qp_; - const VectorXd& dual_solution_; - const VectorXd& dual_gradient_; -}; - // JointTrustRegionProblem defines the joint primal/dual trust region problem // given a QuadraticProgram, primal and dual solutions, primal and dual // gradients, and the primal weight. The joint problem (implicitly) concatenates @@ -242,69 +167,6 @@ struct TrustRegionResultStepSize { double objective_value; }; -// The distance (in the indexed element) from the center point to the bound, in -// the direction that reduces the objective. -template -double DistanceAtCriticalStepSize(const TrustRegionProblem& problem, - const int64_t index) { - if (problem.Objective(index) == 0.0) { - return 0.0; - } - if (problem.Objective(index) > 0.0) { - return problem.LowerBound(index) - problem.CenterPoint(index); - } else { - return problem.UpperBound(index) - problem.CenterPoint(index); - } -} - -// The critical step size is the step size at which the indexed element hits its -// bound (or infinity if that doesn't happen). -template -double CriticalStepSize(const TrustRegionProblem& problem, - const int64_t index) { - if (problem.Objective(index) == 0.0) { - return std::numeric_limits::infinity(); - } - return -problem.NormWeight(index) * - DistanceAtCriticalStepSize(problem, index) / problem.Objective(index); -} - -// The value of the indexed element at the given step_size, projected onto the -// bounds. -template -double ProjectedValue(const TrustRegionProblem& problem, const int64_t index, - const double step_size) { - const double full_step = - problem.CenterPoint(index) - - step_size * problem.Objective(index) / problem.NormWeight(index); - return std::min(std::max(full_step, problem.LowerBound(index)), - problem.UpperBound(index)); -} - -// The distance (in the indexed element) from the center point to the value at -// at the given step size projected onto the bounds. -template -double DistanceAtProjectedValue(const TrustRegionProblem& problem, - const int64_t index, const double step_size) { - return ProjectedValue(problem, index, step_size) - problem.CenterPoint(index); -} - -// This is an easy way of computing medians that's slightly off when the length -// of the array is even. "array" is intentionally passed by copy. -// "value_function" maps an element of "array" to its (double) value. Returns -// the value of the median element. -template -double EasyMedian(ArrayType array, ValueFunction value_function) { - CHECK_GT(array.size(), 0); - auto middle = array.begin() + (array.size() / 2); - absl::c_nth_element(array, middle, - [&](typename ArrayType::value_type lhs, - typename ArrayType::value_type rhs) { - return value_function(lhs) < value_function(rhs); - }); - return value_function(*middle); -} - // "problem" is sharded according to the given sharder. Within each shard, // this function selects the subset of elements corresponding to // indexed_components_by_shard, and takes the median of the critical step sizes @@ -322,9 +184,9 @@ double MedianOfShardMedians( const auto& indexed_shard_components = indexed_components_by_shard[shard.Index()]; if (!indexed_shard_components.empty()) { - shard_medians[shard.Index()] = - EasyMedian(indexed_shard_components, [&](const int64_t index) { - return CriticalStepSize(problem, index); + shard_medians[shard.Index()] = internal::EasyMedian( + indexed_shard_components, [&](const int64_t index) { + return internal::CriticalStepSize(problem, index); }); } }); @@ -335,122 +197,63 @@ double MedianOfShardMedians( } } CHECK(!non_empty_medians.empty()); - return EasyMedian(non_empty_medians, [](const double x) { return x; }); + return internal::EasyMedian(non_empty_medians, + [](const double x) { return x; }); } struct InitialState { std::vector> undecided_components_by_shard; - double radius_coefficient_of_decided_coefficients; + double radius_coefficient_of_decided_components; }; -// Lists the undecided components (per shard) as those with finite critical step -// sizes. The components with infinite critical step sizes will never hit their -// bounds, so we compute their contribution to the radius as -// radius_coefficient_of_decided_coefficients. template InitialState ComputeInitialState(const TrustRegionProblem& problem, const Sharder& sharder) { InitialState result; result.undecided_components_by_shard.resize(sharder.NumShards()); - result.radius_coefficient_of_decided_coefficients = + result.radius_coefficient_of_decided_components = sharder.ParallelSumOverShards([&](const Sharder::Shard& shard) { - std::vector& undecided_components = - result.undecided_components_by_shard[shard.Index()]; const int64_t shard_start = sharder.ShardStart(shard.Index()); const int64_t shard_size = sharder.ShardSize(shard.Index()); - undecided_components.reserve(shard_size); - double radius_coefficient = 0.0; - for (int64_t index = shard_start; index < shard_start + shard_size; - ++index) { - if (std::isfinite(CriticalStepSize(problem, index))) { - undecided_components.push_back(index); - } else { - // Simplified from norm_weight * (objective / norm_weight)^2. - radius_coefficient += MathUtil::Square(problem.Objective(index)) / - problem.NormWeight(index); - } - } - return radius_coefficient; + return internal::ComputeInitialUndecidedComponents( + problem, shard_start, shard_start + shard_size, + result.undecided_components_by_shard[shard.Index()]); }); return result; } template -double RadiusSquaredAtUndecidedComponents( - const TrustRegionProblem& problem, - const std::vector>& undecided_components_by_shard, - const double step_size, const Sharder& sharder) { +double RadiusSquaredOfUndecidedComponents( + const TrustRegionProblem& problem, const double step_size, + const Sharder& sharder, + const std::vector>& undecided_components_by_shard) { return sharder.ParallelSumOverShards([&](const Sharder::Shard& shard) { - const std::vector& undecided_components = - undecided_components_by_shard[shard.Index()]; - return absl::c_accumulate( - undecided_components, 0.0, [&](double sum, int64_t index) { - return sum + problem.NormWeight(index) * - MathUtil::Square(DistanceAtProjectedValue( - problem, index, step_size)); - }); + return internal::RadiusSquaredOfUndecidedComponents( + problem, step_size, undecided_components_by_shard[shard.Index()]); }); } -// Points whose critical step-sizes are greater than or equal to pivot are -// eliminated from the undecided components (we know they'll be determined by -// center_point - step_size * objective / norm_weights). Returns the coefficient -// of step_size^2 that accounts of the contribution of the removed variables -// to the radius squared. template -double RemoveCriticalStepsAbovePivot( - const double pivot, const TrustRegionProblem& problem, +double RemoveCriticalStepsAboveThreshold( + const TrustRegionProblem& problem, const double step_size_threshold, const Sharder& sharder, std::vector>& undecided_components_by_shard) { return sharder.ParallelSumOverShards([&](const Sharder::Shard& shard) { - std::vector& undecided_components = - undecided_components_by_shard[shard.Index()]; - double variable_radius_coefficient = 0.0; - for (const int64_t index : undecided_components) { - if (CriticalStepSize(problem, index) >= pivot) { - // Simplified from norm_weight * (objective / norm_weight)^2. - variable_radius_coefficient += - MathUtil::Square(problem.Objective(index)) / - problem.NormWeight(index); - } - } - auto result = - std::remove_if(undecided_components.begin(), undecided_components.end(), - [&](const int64_t index) { - return CriticalStepSize(problem, index) >= pivot; - }); - undecided_components.erase(result, undecided_components.end()); - return variable_radius_coefficient; + return internal::RemoveCriticalStepsAboveThreshold( + problem, step_size_threshold, + undecided_components_by_shard[shard.Index()]); }); } -// Points whose critical step-sizes are smaller than or equal to pivot are -// eliminated from the undecided components (we know they'll always be at their -// bounds). Returns the weighted distance squared from the center point for the -// removed components. template -double RemoveCriticalStepsBelowPivot( - const double pivot, const TrustRegionProblem& problem, +double RemoveCriticalStepsBelowThreshold( + const TrustRegionProblem& problem, const double step_size_threshold, const Sharder& sharder, std::vector>& undecided_components_by_shard) { return sharder.ParallelSumOverShards([&](const Sharder::Shard& shard) { - std::vector& undecided_components = - undecided_components_by_shard[shard.Index()]; - double radius_sq = 0.0; - for (const int64_t index : undecided_components) { - if (CriticalStepSize(problem, index) <= pivot) { - radius_sq += - problem.NormWeight(index) * - MathUtil::Square(DistanceAtCriticalStepSize(problem, index)); - } - } - auto result = - std::remove_if(undecided_components.begin(), undecided_components.end(), - [&](const int64_t index) { - return CriticalStepSize(problem, index) <= pivot; - }); - undecided_components.erase(result, undecided_components.end()); - return radius_sq; + return internal::RemoveCriticalStepsBelowThreshold( + problem, step_size_threshold, + undecided_components_by_shard[shard.Index()]); }); } @@ -481,7 +284,7 @@ VectorXd ComputeSolution(const TrustRegionProblem& problem, const int64_t shard_size = sharder.ShardSize(shard.Index()); for (int64_t index = shard_start; index < shard_start + shard_size; ++index) { - solution[index] = ProjectedValue(problem, index, step_size); + solution[index] = internal::ProjectedValue(problem, index, step_size); } }); return solution; @@ -497,7 +300,7 @@ double ComputeObjectiveValue(const TrustRegionProblem& problem, for (int64_t index = shard_start; index < shard_start + shard_size; ++index) { shard_value += problem.Objective(index) * - (ProjectedValue(problem, index, step_size) - + (internal::ProjectedValue(problem, index, step_size) - problem.CenterPoint(index)); } return shard_value; @@ -572,7 +375,7 @@ TrustRegionResultStepSize SolveTrustRegionStepSize( // center_point - step_size * objective / norm_weights. These variables are // not at their bounds in the solution, except in degenerate cases. double variable_radius_coefficient = - initial_state.radius_coefficient_of_decided_coefficients; + initial_state.radius_coefficient_of_decided_components; // For each shard, the components of the variables that aren't accounted for // in fixed_radius_squared or variable_radius_coefficient, i.e., we don't know @@ -595,23 +398,23 @@ TrustRegionResultStepSize SolveTrustRegionStepSize( actual_elements_seen += NumUndecidedComponents(undecided_components_by_shard); - const double pivot = + const double step_size_threshold = MedianOfShardMedians(problem, undecided_components_by_shard, sharder); - const double radius_squared_at_undecided_components = - RadiusSquaredAtUndecidedComponents(problem, - undecided_components_by_shard, - /*step_size=*/pivot, sharder); + const double radius_squared_of_undecided_components = + RadiusSquaredOfUndecidedComponents( + problem, /*step_size=*/step_size_threshold, sharder, + undecided_components_by_shard); - const double radius_squared_at_pivot = - radius_squared_at_undecided_components + fixed_radius_squared + - variable_radius_coefficient * MathUtil::Square(pivot); + const double radius_squared_at_threshold = + radius_squared_of_undecided_components + fixed_radius_squared + + variable_radius_coefficient * MathUtil::Square(step_size_threshold); - if (radius_squared_at_pivot > MathUtil::Square(target_radius)) { - variable_radius_coefficient += RemoveCriticalStepsAbovePivot( - pivot, problem, sharder, undecided_components_by_shard); + if (radius_squared_at_threshold > MathUtil::Square(target_radius)) { + variable_radius_coefficient += RemoveCriticalStepsAboveThreshold( + problem, step_size_threshold, sharder, undecided_components_by_shard); } else { - fixed_radius_squared += RemoveCriticalStepsBelowPivot( - pivot, problem, sharder, undecided_components_by_shard); + fixed_radius_squared += RemoveCriticalStepsBelowThreshold( + problem, step_size_threshold, sharder, undecided_components_by_shard); } } VLOG(1) << "Total passes through variables: " @@ -997,8 +800,8 @@ MaxNormBoundResult ComputeMaxNormPrimalTrustRegionBound( const double primal_radius, const VectorXd& dual_product) { LagrangianPart primal_part = ComputePrimalGradient(sharded_qp, primal_solution, dual_product); - PrimalTrustRegionProblem primal_problem(&sharded_qp.Qp(), &primal_solution, - &primal_part.gradient); + internal::PrimalTrustRegionProblem primal_problem( + &sharded_qp.Qp(), &primal_solution, &primal_part.gradient); TrustRegionResultStepSize trust_region_result = SolveTrustRegionStepSize( primal_problem, primal_radius, sharded_qp.PrimalSharder()); return {.part_of_lagrangian_value = primal_part.value, @@ -1010,8 +813,8 @@ MaxNormBoundResult ComputeMaxNormDualTrustRegionBound( const double dual_radius, const VectorXd& primal_product) { LagrangianPart dual_part = ComputeDualGradient(sharded_qp, dual_solution, primal_product); - DualTrustRegionProblem dual_problem(&sharded_qp.Qp(), &dual_solution, - &dual_part.gradient); + internal::DualTrustRegionProblem dual_problem( + &sharded_qp.Qp(), &dual_solution, &dual_part.gradient); TrustRegionResultStepSize trust_region_result = SolveTrustRegionStepSize( dual_problem, dual_radius, sharded_qp.DualSharder()); return {.part_of_lagrangian_value = dual_part.value, diff --git a/ortools/pdlp/trust_region.h b/ortools/pdlp/trust_region.h index 0e07efab24..05b6ab399f 100644 --- a/ortools/pdlp/trust_region.h +++ b/ortools/pdlp/trust_region.h @@ -14,7 +14,16 @@ #ifndef PDLP_TRUST_REGION_H_ #define PDLP_TRUST_REGION_H_ +#include +#include +#include +#include +#include + #include "Eigen/Core" +#include "absl/algorithm/container.h" +#include "ortools/base/logging.h" +#include "ortools/base/mathutil.h" #include "ortools/pdlp/sharded_quadratic_program.h" #include "ortools/pdlp/sharder.h" @@ -162,6 +171,248 @@ LocalizedLagrangianBounds ComputeLocalizedLagrangianBounds( bool use_diagonal_qp_trust_region_solver, double diagonal_qp_trust_region_solver_tolerance); +namespace internal { + +// These functions templated on TrustRegionProblem compute values useful to the +// trust region solve. The templated TrustRegionProblem type should provide +// methods: +// double Objective(int64_t index) const; +// double LowerBound(int64_t index) const; +// double UpperBound(int64_t index) const; +// double CenterPoint(int64_t index) const; +// double NormWeight(int64_t index) const; +// See trust_region.cc for more details and several implementations. + +// The distance (in the indexed element) from the center point to the bound, in +// the direction that reduces the objective. +template +double DistanceAtCriticalStepSize(const TrustRegionProblem& problem, + const int64_t index) { + if (problem.Objective(index) == 0.0) { + return 0.0; + } + if (problem.Objective(index) > 0.0) { + return problem.LowerBound(index) - problem.CenterPoint(index); + } else { + return problem.UpperBound(index) - problem.CenterPoint(index); + } +} + +// The critical step size is the step size at which the indexed element hits its +// bound (or infinity if that doesn't happen). +template +double CriticalStepSize(const TrustRegionProblem& problem, + const int64_t index) { + if (problem.Objective(index) == 0.0) { + return std::numeric_limits::infinity(); + } + return -problem.NormWeight(index) * + DistanceAtCriticalStepSize(problem, index) / problem.Objective(index); +} + +// The value of the indexed element at the given step_size, projected onto the +// bounds. +template +double ProjectedValue(const TrustRegionProblem& problem, const int64_t index, + const double step_size) { + const double full_step = + problem.CenterPoint(index) - + step_size * problem.Objective(index) / problem.NormWeight(index); + return std::min(std::max(full_step, problem.LowerBound(index)), + problem.UpperBound(index)); +} + +// An easy way of computing medians that's slightly off when the length of the +// array is even. "array" is intentionally passed by copy. +// "value_function" maps an element of "array" to its (double) value. Returns +// the value of the median element. +template +double EasyMedian(ArrayType array, ValueFunction value_function) { + CHECK_GT(array.size(), 0); + auto middle = array.begin() + (array.size() / 2); + absl::c_nth_element(array, middle, + [&](typename ArrayType::value_type lhs, + typename ArrayType::value_type rhs) { + return value_function(lhs) < value_function(rhs); + }); + return value_function(*middle); +} + +// Lists the undecided components (from [start_index, end_index) as those with +// finite critical step sizes. The components with infinite critical step sizes +// will never hit their bounds, so returns their contribution to square of the +// radius. +template +double ComputeInitialUndecidedComponents( + const TrustRegionProblem& problem, int64_t start_index, int64_t end_index, + std::vector& undecided_components) { + // TODO(user): Evaluate dropping this reserve(), since it wastes space + // if many components are decided. + undecided_components.clear(); + undecided_components.reserve(end_index - start_index); + double radius_coefficient = 0.0; + for (int64_t index = start_index; index < end_index; ++index) { + if (std::isfinite(internal::CriticalStepSize(problem, index))) { + undecided_components.push_back(index); + } else { + // Simplified from norm_weight * (objective / norm_weight)^2. + radius_coefficient += MathUtil::Square(problem.Objective(index)) / + problem.NormWeight(index); + } + } + return radius_coefficient; +} + +template +double RadiusSquaredOfUndecidedComponents( + const TrustRegionProblem& problem, const double step_size, + const std::vector& undecided_components) { + return absl::c_accumulate( + undecided_components, 0.0, [&](double sum, int64_t index) { + const double distance_at_projected_value = + internal::ProjectedValue(problem, index, step_size) - + problem.CenterPoint(index); + return sum + problem.NormWeight(index) * + MathUtil::Square(distance_at_projected_value); + }); +} + +// Points whose critical step-sizes are greater than or equal to +// step_size_threshold are eliminated from the undecided components (we know +// they'll be determined by center_point - step_size * objective / +// norm_weights). Returns the coefficient of step_size^2 that accounts of the +// contribution of the removed variables to the radius squared. +template +double RemoveCriticalStepsAboveThreshold( + const TrustRegionProblem& problem, const double step_size_threshold, + std::vector& undecided_components) { + double variable_radius_coefficient = 0.0; + for (const int64_t index : undecided_components) { + if (internal::CriticalStepSize(problem, index) >= step_size_threshold) { + // Simplified from norm_weight * (objective / norm_weight)^2. + variable_radius_coefficient += + MathUtil::Square(problem.Objective(index)) / + problem.NormWeight(index); + } + } + auto result = + std::remove_if(undecided_components.begin(), undecided_components.end(), + [&](const int64_t index) { + return internal::CriticalStepSize(problem, index) >= + step_size_threshold; + }); + undecided_components.erase(result, undecided_components.end()); + return variable_radius_coefficient; +} + +// Points whose critical step-sizes are smaller than or equal to +// step_size_threshold are eliminated from the undecided components (we know +// they'll always be at their bounds). Returns the weighted distance squared +// from the center point for the removed components. +template +double RemoveCriticalStepsBelowThreshold( + const TrustRegionProblem& problem, const double step_size_threshold, + std::vector& undecided_components) { + double radius_sq = 0.0; + for (const int64_t index : undecided_components) { + if (internal::CriticalStepSize(problem, index) <= step_size_threshold) { + radius_sq += problem.NormWeight(index) * + MathUtil::Square( + internal::DistanceAtCriticalStepSize(problem, index)); + } + } + auto result = + std::remove_if(undecided_components.begin(), undecided_components.end(), + [&](const int64_t index) { + return internal::CriticalStepSize(problem, index) <= + step_size_threshold; + }); + undecided_components.erase(result, undecided_components.end()); + return radius_sq; +} + +// PrimalTrustRegionProblem defines the primal trust region problem given a +// QuadraticProgram, primal solution, and primal gradient. It captures const +// references to the constructor arguments, which should outlive the class +// instance. +// The corresponding trust region problem is +// min_x primal_gradient' * (x - primal_solution) +// s.t. qp.variable_lower_bounds <= x <= qp.variable_upper_bounds +// || x - primal_solution ||_2 <= target_radius +class PrimalTrustRegionProblem { + public: + PrimalTrustRegionProblem(const QuadraticProgram* qp, + const Eigen::VectorXd* primal_solution, + const Eigen::VectorXd* primal_gradient, + const double norm_weight = 1.0) + : qp_(*qp), + primal_solution_(*primal_solution), + primal_gradient_(*primal_gradient), + norm_weight_(norm_weight) {} + double Objective(int64_t index) const { return primal_gradient_[index]; } + double LowerBound(int64_t index) const { + return qp_.variable_lower_bounds[index]; + } + double UpperBound(int64_t index) const { + return qp_.variable_upper_bounds[index]; + } + double CenterPoint(int64_t index) const { return primal_solution_[index]; } + double NormWeight(int64_t index) const { return norm_weight_; } + + private: + const QuadraticProgram& qp_; + const Eigen::VectorXd& primal_solution_; + const Eigen::VectorXd& primal_gradient_; + const double norm_weight_; +}; + +// DualTrustRegionProblem defines the dual trust region problem given a +// QuadraticProgram, dual solution, and dual gradient. It captures const +// references to the constructor arguments, which should outlive the class +// instance. +// The corresponding trust region problem is +// max_y dual_gradient' * (y - dual_solution) +// s.t. qp.implicit_dual_lower_bounds <= y <= qp.implicit_dual_upper_bounds +// || y - dual_solution ||_2 <= target_radius +// where the implicit dual bounds are those given in +// https://developers.google.com/optimization/lp/pdlp_math#dual_variable_bounds +class DualTrustRegionProblem { + public: + DualTrustRegionProblem(const QuadraticProgram* qp, + const Eigen::VectorXd* dual_solution, + const Eigen::VectorXd* dual_gradient, + const double norm_weight = 1.0) + : qp_(*qp), + dual_solution_(*dual_solution), + dual_gradient_(*dual_gradient), + norm_weight_(norm_weight) {} + double Objective(int64_t index) const { + // The objective is negated because the trust region problem objective is + // minimize, but for the dual problem we want to maximize the gradient. + return -dual_gradient_[index]; + } + double LowerBound(int64_t index) const { + return std::isfinite(qp_.constraint_upper_bounds[index]) + ? -std::numeric_limits::infinity() + : 0.0; + } + double UpperBound(int64_t index) const { + return std::isfinite(qp_.constraint_lower_bounds[index]) + ? std::numeric_limits::infinity() + : 0.0; + } + double CenterPoint(int64_t index) const { return dual_solution_[index]; } + double NormWeight(int64_t index) const { return norm_weight_; } + + private: + const QuadraticProgram& qp_; + const Eigen::VectorXd& dual_solution_; + const Eigen::VectorXd& dual_gradient_; + const double norm_weight_; +}; + +} // namespace internal + } // namespace operations_research::pdlp #endif // PDLP_TRUST_REGION_H_ diff --git a/ortools/pdlp/trust_region_test.cc b/ortools/pdlp/trust_region_test.cc index 5fe5d3c0aa..3634c76e86 100644 --- a/ortools/pdlp/trust_region_test.cc +++ b/ortools/pdlp/trust_region_test.cc @@ -646,19 +646,23 @@ TEST_P(ComputeLocalizedLagrangianBoundsTest, OptimalInBoundRange) { primal_solution << 0.0, 0.0, 0.0, 3.0; VectorXd dual_solution = 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 / 6.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=*/nullptr, /*dual_product=*/nullptr, use_diagonal_qp_solver, /*diagonal_qp_trust_region_solver_tolerance=*/1.0e-6); @@ -690,30 +694,40 @@ TEST_P(ComputeLocalizedLagrangianBoundsTest, OptimalNotInBoundRange) { const double expected_lagrangian = 3.0; EXPECT_DOUBLE_EQ(bounds.lagrangian_value, expected_lagrangian); - // The rough intervals below are based on the observations that: - // 1) All components of the gradients are between 1 and 10. - // 2) The radius of 0.1 is small enough that the variable bounds have a minor - // effect. - // 3) For the trust-region problem with no variable bounds, we have that: - // objective_vector' * (x - center_point) = -target_radius * - // ||objective_vector||. + // Because the dual solution is all zero, the primal gradient is just the + // objective, [5.5, -2, -1, 1]. The dual gradient is the dual subgradient + // coefficient minus the primal product. With a zero dual, for one-sided + // constraints, the dual subgradient coefficient is the bound, and for + // two-sided constraints it is the violated bound (or zero if feasible). Thus, + // the dual subgradient coefficients are [12, 7, -4, -1], and the primal + // product is [6, 0, 0, -3], giving a dual gradient of [6, 7, -4, 2]. switch (primal_dual_norm) { case PrimalDualNorm::kMaxNorm: - // ||objective_vector|| is between 2 and 20 (because the vector has - // length 4), and target_radius = sqrt(2) * 0.1 ≈ 0.14. - // So the bounds should move between 0.28 and 2.8. - EXPECT_LE(bounds.lower_bound, expected_lagrangian - 0.28); - EXPECT_GE(bounds.lower_bound, expected_lagrangian - 2.8); - EXPECT_GE(bounds.upper_bound, expected_lagrangian + 0.28); - EXPECT_LE(bounds.upper_bound, expected_lagrangian + 2.8); + // The target_radius r = sqrt(2) * 0.1 ≈ 0.14, and the projected primal + // direction is d=[-5.5, 2, 1, -1]. The resulting delta is d / ||d|| * r, + // giving an objective delta of + // ||d|| * r. + EXPECT_NEAR(bounds.lower_bound, + expected_lagrangian - 0.1 * sqrt(2) * sqrt(36.25), 1.0e-6); + // The target_radius r = sqrt(2) * 0.1 ≈ 0.14, and the projected dual + // direction is d=[6, 0, 0, 2]. The resulting delta is d / ||d|| * r, + // giving an objective delta of + // ||d|| * r. + EXPECT_NEAR(bounds.upper_bound, + expected_lagrangian + 0.1 * sqrt(2) * sqrt(40.0), 1.0e-6); break; case PrimalDualNorm::kEuclideanNorm: - // In this case, the value of objective_vector' * (x - center_point) is - // upper_bound - lower_bound. ||objective_vector|| is between 2.8 and 28, - // and target_radius ≈ 0.14. - EXPECT_GE(bounds.upper_bound - bounds.lower_bound, 0.39); - EXPECT_LE(bounds.upper_bound - bounds.lower_bound, 3.9); + // In this case, r = target_radius * sqrt(2) (because the euclidean norm + // includes a factor of 0.5). The projected combined direction is d=[-5.5, + // 2, 1, -1; 6, 0, 0, 2]. The resulting primal delta is d[primal] / ||d|| + // * r, and the resulting dual delta is d[dual] / ||d|| * r. + EXPECT_NEAR(bounds.lower_bound, + expected_lagrangian - 0.1 * sqrt(2) * 36.25 / sqrt(76.25), + 1.0e-6); + EXPECT_NEAR(bounds.upper_bound, + expected_lagrangian + 0.1 * sqrt(2) * 40 / sqrt(76.25), + 1.0e-6); break; } } diff --git a/ortools/sat/encoding.cc b/ortools/sat/encoding.cc index d54628980a..174ac1556a 100644 --- a/ortools/sat/encoding.cc +++ b/ortools/sat/encoding.cc @@ -31,10 +31,7 @@ namespace operations_research { namespace sat { EncodingNode::EncodingNode(Literal l) - : depth_(0), - lb_(0), - ub_(1), - for_sorting_(l.Variable()), + : for_sorting_(l.Variable()), child_a_(nullptr), child_b_(nullptr), literals_(1, l) {} @@ -75,6 +72,21 @@ void EncodingNode::InitializeLazyNode(EncodingNode* a, EncodingNode* b, for_sorting_ = std::min(a->for_sorting_, b->for_sorting_); } +void EncodingNode::InitializeLazyCoreNode(Coefficient weight, EncodingNode* a, + EncodingNode* b) { + CHECK(literals_.empty()) << "Already initialized"; + child_a_ = a; + child_b_ = b; + ub_ = a->ub_ + b->ub_; + weight_ = weight; + weight_lb_ = a->lb_ + b->lb_; + lb_ = weight_lb_ + 1; + depth_ = 1 + std::max(a->depth_, b->depth_); + + // Merging the node of the same depth in order seems to help a bit. + for_sorting_ = std::min(a->for_sorting_, b->for_sorting_); +} + bool EncodingNode::IncreaseCurrentUB(SatSolver* solver) { if (current_ub() == ub_) return false; literals_.emplace_back(BooleanVariable(solver->NumVariables()), true); @@ -86,7 +98,7 @@ bool EncodingNode::IncreaseCurrentUB(SatSolver* solver) { return true; } -int EncodingNode::Reduce(const SatSolver& solver) { +Coefficient EncodingNode::Reduce(const SatSolver& solver) { int i = 0; while (i < literals_.size() && solver.Assignment().LiteralIsTrue(literals_[i])) { @@ -99,24 +111,54 @@ int EncodingNode::Reduce(const SatSolver& solver) { literals_.pop_back(); ub_ = lb_ + literals_.size(); } - return i; + + if (weight_lb_ >= lb_) return Coefficient(0); + const Coefficient result = Coefficient(lb_ - weight_lb_) * weight_; + weight_lb_ = lb_; + return result; } -void EncodingNode::ApplyUpperBound(int64_t upper_bound, SatSolver* solver) { - if (size() <= upper_bound) return; - for (int i = upper_bound; i < size(); ++i) { +void EncodingNode::ApplyWeightUpperBound(Coefficient gap, SatSolver* solver) { + CHECK_GT(weight_, 0); + const Coefficient num_allowed = (gap / weight_); + const Coefficient new_size = + std::max(Coefficient(0), Coefficient(weight_lb_ - lb_) + num_allowed); + if (size() <= new_size) return; + for (int i = new_size.value(); i < size(); ++i) { solver->AddUnitClause(literal(i).Negated()); } - literals_.resize(upper_bound); + literals_.resize(new_size.value()); ub_ = lb_ + literals_.size(); } +Literal EncodingNode::GetAssumption(SatSolver* solver) { + CHECK(!HasNoWeight()); + const int index = weight_lb_ - lb_; + CHECK_GE(index, 0) << "Not reduced?"; + while (index >= literals_.size()) { + CHECK_NE(solver, nullptr); + IncreaseNodeSize(this, solver); + } + return literals_[index].Negated(); +} + +void EncodingNode::IncreaseWeightLb() { + CHECK_LT(weight_lb_ - lb_, literals_.size()); + weight_lb_++; +} + +bool EncodingNode::HasNoWeight() const { + return weight_ == 0 || weight_lb_ >= ub_; +} + std::string EncodingNode::DebugString( const VariablesAssignment& assignment) const { std::string result; absl::StrAppend(&result, "depth:", depth_); absl::StrAppend(&result, " [", lb_, ",", lb_ + literals_.size(), "]"); absl::StrAppend(&result, " ub:", ub_); + absl::StrAppend(&result, " weight:", weight_.value()); + absl::StrAppend(&result, " weight_lb:", weight_lb_); absl::StrAppend(&result, " values:"); const size_t limit = 20; int value = 0; @@ -312,13 +354,13 @@ struct SortEncodingNodePointers { }; } // namespace -EncodingNode* LazyMergeAllNodeWithPQ(const std::vector& nodes, - SatSolver* solver, - std::deque* repository) { +EncodingNode* LazyMergeAllNodeWithPQAndIncreaseLb( + Coefficient weight, const std::vector& nodes, + SatSolver* solver, std::deque* repository) { std::priority_queue, SortEncodingNodePointers> pq(nodes.begin(), nodes.end()); - while (pq.size() > 1) { + while (pq.size() > 2) { EncodingNode* a = pq.top(); pq.pop(); EncodingNode* b = pq.top(); @@ -326,7 +368,18 @@ EncodingNode* LazyMergeAllNodeWithPQ(const std::vector& nodes, repository->push_back(LazyMerge(a, b, solver)); pq.push(&repository->back()); } - return pq.top(); + + CHECK_EQ(pq.size(), 2); + EncodingNode* a = pq.top(); + pq.pop(); + EncodingNode* b = pq.top(); + pq.pop(); + + repository->push_back(EncodingNode()); + EncodingNode* n = &repository->back(); + n->InitializeLazyCoreNode(weight, a, b); + solver->AddBinaryClause(a->literal(0), b->literal(0)); + return n; } std::vector CreateInitialEncodingNodes( @@ -389,8 +442,6 @@ bool EncodingNodeByDepth(const EncodingNode* a, const EncodingNode* b) { return a->depth() < b->depth(); } -bool EmptyEncodingNode(const EncodingNode* a) { return a->size() == 0; } - } // namespace std::vector ReduceNodesAndExtractAssumptions( @@ -402,14 +453,7 @@ std::vector ReduceNodesAndExtractAssumptions( // at the root node in order to work. solver->Backtrack(0); for (EncodingNode* n : *nodes) { - *lower_bound += n->Reduce(*solver) * n->weight(); - - // Tricky: When we solve a partial set of assumptions, some unused nodes - // might have literals that get fixed to one. If that happens we do need - // to create new literals so we can fix these nodes to their lower bound. - if (n->size() == 0 && n->lb() < n->ub()) { - IncreaseNodeSize(n, solver); - } + *lower_bound += n->Reduce(*solver); } // Fix the nodes right-most variables that are above the gap. @@ -418,12 +462,13 @@ std::vector ReduceNodesAndExtractAssumptions( const Coefficient gap = upper_bound - *lower_bound; if (gap < 0) return {}; for (EncodingNode* n : *nodes) { - n->ApplyUpperBound((gap / n->weight()).value(), solver); + n->ApplyWeightUpperBound(gap, solver); } } // Remove the empty nodes. - nodes->erase(std::remove_if(nodes->begin(), nodes->end(), EmptyEncodingNode), + nodes->erase(std::remove_if(nodes->begin(), nodes->end(), + [](EncodingNode* a) { return a->HasNoWeight(); }), nodes->end()); // Sort the nodes. @@ -447,7 +492,7 @@ std::vector ReduceNodesAndExtractAssumptions( std::vector assumptions; for (EncodingNode* n : *nodes) { if (n->weight() >= stratified_lower_bound) { - assumptions.push_back(n->literal(0).Negated()); + assumptions.push_back(n->GetAssumption(solver)); } } return assumptions; @@ -458,8 +503,8 @@ Coefficient ComputeCoreMinWeight(const std::vector& nodes, Coefficient min_weight = kCoefficientMax; int index = 0; for (int i = 0; i < core.size(); ++i) { - for (; - index < nodes.size() && nodes[index]->literal(0).Negated() != core[i]; + for (; index < nodes.size() && + nodes[index]->GetAssumption(nullptr) != core[i]; ++index) { } CHECK_LT(index, nodes.size()); @@ -485,19 +530,8 @@ bool ProcessCore(const std::vector& core, Coefficient min_weight, std::vector* nodes, SatSolver* solver) { // Backtrack to be able to add new constraints. solver->ResetToLevelZero(); - if (core.size() == 1) { - if (!solver->AddUnitClause(core[0].Negated())) return false; - - // The core will be reduced at the beginning of the next loop. - // Find the associated node, and call IncreaseNodeSize() on it. - for (EncodingNode* n : *nodes) { - if (n->literal(0).Negated() == core[0]) { - IncreaseNodeSize(n, solver); - return true; - } - } - LOG(FATAL) << "Node with literal " << core[0] << " not found!"; + return solver->AddUnitClause(core[0].Negated()); } // Remove from nodes the EncodingNode in the core, merge them, and add the @@ -509,7 +543,7 @@ bool ProcessCore(const std::vector& core, Coefficient min_weight, // Since the nodes appear in order in the core, we can find the // relevant "objective" variable efficiently with a simple linear scan // in the nodes vector (done with index). - for (; (*nodes)[index]->literal(0).Negated() != core[i]; ++index) { + for (; (*nodes)[index]->GetAssumption(nullptr) != core[i]; ++index) { CHECK_LT(index, nodes->size()); (*nodes)[new_node_index] = (*nodes)[index]; ++new_node_index; @@ -534,10 +568,78 @@ bool ProcessCore(const std::vector& core, Coefficient min_weight, ++new_node_index; } nodes->resize(new_node_index); - nodes->push_back(LazyMergeAllNodeWithPQ(to_merge, solver, repository)); - IncreaseNodeSize(nodes->back(), solver); - nodes->back()->set_weight(min_weight); - return solver->AddUnitClause(nodes->back()->literal(0)); + nodes->push_back(LazyMergeAllNodeWithPQAndIncreaseLb(min_weight, to_merge, + solver, repository)); + return !solver->IsModelUnsat(); +} + +bool ProcessCoreWithNewEncoding(const std::vector& core, + Coefficient min_weight, + std::deque* repository, + std::vector* nodes, + SatSolver* solver) { + // Backtrack to be able to add new constraints. + solver->ResetToLevelZero(); + + if (core.size() == 1) { + return solver->AddUnitClause(core[0].Negated()); + } + + std::vector new_nodes; + std::vector to_merge; + + // Preconditions. + for (EncodingNode* n : *nodes) { + CHECK_GT(n->size(), 0); + } + + // Remove from nodes the EncodingNode in the core, merge them, and add the + // resulting EncodingNode at the back. + int index = 0; + for (int i = 0; i < core.size(); ++i) { + // Since the nodes appear in order in the core, we can find the + // relevant "objective" variable efficiently with a simple linear scan + // in the nodes vector (done with index). + CHECK_LT(index, nodes->size()); + for (; (*nodes)[index]->GetAssumption(nullptr) != core[i]; ++index) { + CHECK_LT(index, nodes->size()); + new_nodes.push_back((*nodes)[index]); + } + CHECK_LT(index, nodes->size()); + + // We have a node from the core. + // We will distinguish its first literal. + EncodingNode* n = (*nodes)[index]; + const Literal lit = n->GetAssumption(nullptr).Negated(); + n->IncreaseWeightLb(); + ++index; + CHECK_GT(n->size(), 0); + + // TODO(user): For node with same depth, the sorting order is not the same + // if we create a new node or reuse one. Experiment what is the best order. + repository->emplace_back(lit); + EncodingNode* new_bool_node = &repository->back(); + new_bool_node->set_depth(n->depth()); + CHECK_GT(new_bool_node->size(), 0); + to_merge.push_back(new_bool_node); + if (n->weight() > min_weight) { + new_bool_node->set_weight(n->weight() - min_weight); + new_nodes.push_back(new_bool_node); + } + + if (!n->HasNoWeight()) { + new_nodes.push_back(n); + } + } + + for (; index < nodes->size(); ++index) { + new_nodes.push_back((*nodes)[index]); + } + + new_nodes.push_back(LazyMergeAllNodeWithPQAndIncreaseLb(min_weight, to_merge, + solver, repository)); + *nodes = new_nodes; + return !solver->IsModelUnsat(); } } // namespace sat diff --git a/ortools/sat/encoding.h b/ortools/sat/encoding.h index 5bba3ed62f..8a8483b54b 100644 --- a/ortools/sat/encoding.h +++ b/ortools/sat/encoding.h @@ -14,6 +14,8 @@ // Algorithms to encode constraints into their SAT representation. Currently, // this contains one possible encoding of a cardinality constraint as used by // the core-based optimization algorithm in optimization.h. +// +// This is also known as the incremental totalizer encoding in the literature. #ifndef OR_TOOLS_SAT_ENCODING_H_ #define OR_TOOLS_SAT_ENCODING_H_ @@ -70,6 +72,8 @@ class EncodingNode { // Only one literals will be created by this operation. Note that no clauses // linking it with a or b are added by this function. void InitializeLazyNode(EncodingNode* a, EncodingNode* b, SatSolver* solver); + void InitializeLazyCoreNode(Coefficient weight, EncodingNode* a, + EncodingNode* b); // Returns a literal with the meaning 'this node number is > i'. // The given i must be in [lb_, current_ub). @@ -94,20 +98,33 @@ class EncodingNode { // Returns false if we were already at the upper bound for this node. bool IncreaseCurrentUB(SatSolver* solver); - // Removes the left-side literals fixed to 1 and returns the number of - // literals removed this way. Note that this increases lb_ and reduces the - // number of active literals. It also removes any right-side literals fixed to - // 0. If such a literal exists, ub is updated accordingly. - int Reduce(const SatSolver& solver); + // Removes the left-side literals fixed to 1. Note that this increases lb_ and + // reduces the number of active literals. It also removes any right-side + // literals fixed to 0. If such a literal exists, ub is updated accordingly. + // + // Return the overall weight increase. + Coefficient Reduce(const SatSolver& solver); - // Fix the right-side variables with indices >= to the given upper_bound to - // false. - void ApplyUpperBound(int64_t upper_bound, SatSolver* solver); + // GetAssumption() might need to create new literals. + Literal GetAssumption(SatSolver* solver); + bool HasNoWeight() const; + void IncreaseWeightLb(); - void set_weight(Coefficient w) { weight_ = w; } + // Fix any literal that would cause the weight of this node to go over the + // gap. + void ApplyWeightUpperBound(Coefficient gap, SatSolver* solver); + + void set_weight(Coefficient w) { + weight_lb_ = lb_; + weight_ = w; + } Coefficient weight() const { return weight_; } + // The depth is mainly used as an heuristic to decide which nodes to merge + // first. See the < operator. + void set_depth(int depth) { depth_ = depth; } int depth() const { return depth_; } + int lb() const { return lb_; } int current_ub() const { return lb_ + literals_.size(); } int ub() const { return ub_; } @@ -118,11 +135,14 @@ class EncodingNode { std::string DebugString(const VariablesAssignment& assignment) const; private: - int depth_; - int lb_; - int ub_; + int depth_ = 0; + int lb_ = 0; + int ub_ = 1; BooleanVariable for_sorting_; + // The weight is only applies for literal >= this lb. + int weight_lb_ = 0; + Coefficient weight_; EncodingNode* child_a_; EncodingNode* child_b_; @@ -131,24 +151,6 @@ class EncodingNode { std::vector literals_; }; -// Note that we use <= because on 32 bits architecture, the size will actually -// be smaller than 64 bytes. One exception is with visual studio on windows, in -// debug mode, where the struct is bigger. -#if defined(_M_X64) && defined(_DEBUG) -// In debug, with msvc, std::Vector is 32 -static_assert(sizeof(EncodingNode) == 72, - "ERROR_EncodingNode_is_not_well_compacted"); -#elif defined(_GLIBCXX_DEBUG) -// In debug GNU libstdc++ std::Vector is up to 56 -static_assert(sizeof(EncodingNode) <= 96, - "ERROR_EncodingNode_is_not_well_compacted"); -#else -// Note that we use <= because on 32 bits architecture, the size will actually -// be smaller than 64 bytes. -static_assert(sizeof(EncodingNode) <= 64, - "ERROR_EncodingNode_is_not_well_compacted"); -#endif - // Merges the two given EncodingNodes by creating a new node that corresponds to // the sum of the two given ones. Only the left-most binary variable is created // for the parent node, the other ones will be created later when needed. @@ -174,10 +176,11 @@ EncodingNode* MergeAllNodesWithDeque(Coefficient upper_bound, std::deque* repository); // Same as MergeAllNodesWithDeque() but use a priority queue to merge in -// priority nodes with smaller sizes. -EncodingNode* LazyMergeAllNodeWithPQ(const std::vector& nodes, - SatSolver* solver, - std::deque* repository); +// priority nodes with smaller sizes. This also enforce that the sum of nodes +// is greater than its lower bound. +EncodingNode* LazyMergeAllNodeWithPQAndIncreaseLb( + Coefficient weight, const std::vector& nodes, + SatSolver* solver, std::deque* repository); // Returns a vector with one new EncodingNode by variable in the given // objective. Sets the offset to the negated sum of the negative coefficient, @@ -216,6 +219,24 @@ bool ProcessCore(const std::vector& core, Coefficient min_weight, std::deque* repository, std::vector* nodes, SatSolver* solver); +// There is more than one way to create new assumptions and encode the +// information from this core. This is slightly different from ProcessCore() and +// follow the algorithm used by many of the top max-SAT solver under the name +// incremental OLL. This is described in: +// António Morgado, Carmine Dodaro, Joao Marques-Silva. "Core-Guided MaxSAT +// with Soft Cardinality Constraints". CP 2014. pp. 564-573. +// António Morgado, Alexey Ignatiev, Joao Marques-Silva. "MSCG: Robust +// Core-Guided MaxSAT Solving." JSAT 9. 2014. pp. 129-134. +// +// TODO(user): The last time this was tested, it was however not as good as the +// ProcessCore() version. That might change as we code/change more heuristic, so +// we keep it around. +bool ProcessCoreWithAlternativeEncoding(const std::vector& core, + Coefficient min_weight, + std::deque* repository, + std::vector* nodes, + SatSolver* solver); + } // namespace sat } // namespace operations_research diff --git a/ortools/sat/optimization.cc b/ortools/sat/optimization.cc index dfb3aef855..1b38ad0066 100644 --- a/ortools/sat/optimization.cc +++ b/ortools/sat/optimization.cc @@ -1625,7 +1625,8 @@ SatSolver::Status CoreBasedOptimizer::OptimizeWithSatEncoding( // The lower bound will be increased by that much. const Coefficient min_weight = ComputeCoreMinWeight(nodes, core); previous_core_info = - absl::StrFormat("core:%u mw:%d", core.size(), min_weight.value()); + absl::StrFormat("core:%u mw:%d d:%d", core.size(), min_weight.value(), + nodes.back()->depth()); // We only count an iter when we found a core. ++iter; diff --git a/ortools/sat/python/swig_helper.cc b/ortools/sat/python/swig_helper.cc new file mode 100644 index 0000000000..3d39b0fe5a --- /dev/null +++ b/ortools/sat/python/swig_helper.cc @@ -0,0 +1,96 @@ +// Copyright 2010-2021 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/sat/swig_helper.h" + +#include + +#include "pybind11/include/pybind11/functional.h" +#include "pybind11/include/pybind11/pybind11.h" +#include "pybind11/include/pybind11/stl.h" +#include "pybind11_protobuf/wrapped_proto_caster.h" +#include "ortools/sat/cp_model.pb.h" +#include "ortools/util/sorted_interval_list.h" + +using ::operations_research::sat::CpSatHelper; +using ::operations_research::sat::CpSolverResponse; +using ::operations_research::sat::SolutionCallback; +using ::operations_research::sat::SolveWrapper; +using ::pybind11::arg; +using ::pybind11_protobuf::WithWrappedProtos; + +class PySolutionCallback : public SolutionCallback { + public: + using SolutionCallback::SolutionCallback; /* Inherit constructors */ + + void OnSolutionCallback() const override { + PYBIND11_OVERRIDE_PURE( + void, /* Return type */ + SolutionCallback, /* Parent class */ + OnSolutionCallback, /* Name of function */ + /* This function has no arguments. The trailing comma + in the previous line is needed for some compilers */ + ); + } +}; + +PYBIND11_MODULE(swig_helper, m) { + pybind11_protobuf::ImportWrappedProtoCasters(); + pybind11::module::import( + + pybind11::class_(m, "SolutionCallback") + .def(pybind11::init<>()) + .def("OnSolutionCallback", &SolutionCallback::OnSolutionCallback) + .def("BestObjectiveBound", &SolutionCallback::BestObjectiveBound) + .def("DeterministicTime", &SolutionCallback::DeterministicTime) + .def("HasResponse", &SolutionCallback::HasResponse) + .def("NumBinaryPropagations", &SolutionCallback::NumBinaryPropagations) + .def("NumBooleans", &SolutionCallback::NumBooleans) + .def("NumBranches", &SolutionCallback::NumBranches) + .def("NumConflicts", &SolutionCallback::NumConflicts) + .def("NumIntegerPropagations", &SolutionCallback::NumIntegerPropagations) + .def("ObjectiveValue", &SolutionCallback::ObjectiveValue) + .def("Response", WithWrappedProtos(&SolutionCallback::Response)) + .def("SolutionBooleanValue", &SolutionCallback::SolutionBooleanValue, + arg("index")) + .def("SolutionIntegerValue", &SolutionCallback::SolutionIntegerValue, + arg("index")) + .def("StopSearch", &SolutionCallback::StopSearch) + .def("UserTime", &SolutionCallback::UserTime) + .def("WallTime", &SolutionCallback::WallTime); + + pybind11::class_(m, "SolveWrapper") + .def(pybind11::init<>()) + .def("AddLogCallback", &SolveWrapper::AddLogCallback, arg("log_callback")) + .def("AddSolutionCallback", &SolveWrapper::AddSolutionCallback, + arg("callback")) + .def("ClearSolutionCallback", &SolveWrapper::ClearSolutionCallback) + .def("SetParameters", WithWrappedProtos(&SolveWrapper::SetParameters), + arg("parameters")) + .def("Solve", WithWrappedProtos(&SolveWrapper::Solve), arg("model_proto")) + .def("StopSearch", &SolveWrapper::StopSearch); + + pybind11::class_(m, "CpSatHelper") + .def_static("ModelStats", WithWrappedProtos(&CpSatHelper::ModelStats)) + .def_static("SolverResponseStats", + WithWrappedProtos(&CpSatHelper::SolverResponseStats)) + .def_static("ValidateModel", + WithWrappedProtos(&CpSatHelper::ValidateModel), + arg("model_proto")) + .def_static("VariableDomain", + WithWrappedProtos(&CpSatHelper::VariableDomain), + arg("variable_proto")) + .def_static("WriteModelToFile", + WithWrappedProtos(&CpSatHelper::WriteModelToFile), + arg("model_proto"), arg("filename")); +} diff --git a/ortools/sat/sat_solver.cc b/ortools/sat/sat_solver.cc index 2deaa5fc6d..37aae1cbc0 100644 --- a/ortools/sat/sat_solver.cc +++ b/ortools/sat/sat_solver.cc @@ -660,7 +660,7 @@ bool SatSolver::PropagateAndStopAfterOneConflictResolution() { // A conflict occurred, compute a nice reason for this failure. same_reason_identifier_.Clear(); const int max_trail_index = ComputeMaxTrailIndex(trail_->FailingClause()); - if (!assumptions_.empty()) { + if (!assumptions_.empty() && !trail_->FailingClause().empty()) { // If the failing clause only contains literal at the assumptions level, // we cannot use the ComputeFirstUIPConflict() code as we might have more // than one decision. @@ -1703,6 +1703,10 @@ void SatSolver::ProcessNewlyFixedVariables() { } // We also clean the binary implication graph. + // Tricky: If we added the first binary clauses above, the binary graph + // is not in "propagated" state as it should be, so we call Propagate() so + // all the checks are happy. + CHECK(binary_implication_graph_->Propagate(trail_)); binary_implication_graph_->RemoveFixedVariables(); num_processed_fixed_variables_ = trail_->Index(); deterministic_time_of_last_fixed_variables_cleanup_ = deterministic_time(); diff --git a/ortools/util/python/sorted_interval_list.cc b/ortools/util/python/sorted_interval_list.cc new file mode 100644 index 0000000000..a676a4e706 --- /dev/null +++ b/ortools/util/python/sorted_interval_list.cc @@ -0,0 +1,45 @@ +// Copyright 2010-2021 Google LLC +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/util/sorted_interval_list.h" + +#include + +#include "pybind11/include/pybind11/pybind11.h" +#include "pybind11/include/pybind11/stl.h" + +using ::operations_research::Domain; +using ::pybind11::arg; + +PYBIND11_MODULE(sorted_interval_list, m) { + pybind11::class_(m, "Domain") + .def_static("AllValues", &Domain::AllValues) + .def_static("FromValues", &Domain::FromValues, arg("values")) + .def_static("FromIntervals", &Domain::FromVectorIntervals, + arg("intervals")) + .def_static("FromFlatIntervals", &Domain::FromFlatIntervals, + arg("flat_intervals")) + .def(pybind11::init()) + .def("AdditionWith", &Domain::AdditionWith, arg("domain")) + .def("Complement", &Domain::Complement) + .def("Contains", &Domain::Contains, arg("value")) + .def("FlattenedIntervals", &Domain::FlattenedIntervals) + .def("IntersectionWith", &Domain::IntersectionWith, arg("domain")) + .def("IsEmpty", &Domain::IsEmpty) + .def("Size", &Domain::Size) + .def("Max", &Domain::Max) + .def("Min", &Domain::Min) + .def("Negation", &Domain::Negation) + .def("UnionWith", &Domain::UnionWith, arg("domain")) + .def("__str__", &Domain::ToString); +}