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

bump pdlp; add log to linear_solver proto; fix bug in CP-SAT

Browse Source
This commit is contained in:
Laurent Perron
2022-02-21 17:26:34 +01:00
parent 8524e60172
commit 951e8bbc1e
21 changed files with 807 additions and 468 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
makefiles/Makefile.dotnet.mk
View File

@@ -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
#############

1
makefiles/Makefile.java.mk
View File

@@ -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)
#############

3
makefiles/Makefile.python.mk
View File

@@ -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*
#############

6
ortools/linear_solver/linear_solver.proto
View File

@@ -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.

49
ortools/linear_solver/pdlp_interface.cc
View File

@@ -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<pdlp::SolutionResponseAndLog> solution =
absl::StatusOr<MPSolutionResponse> 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<MPSolver::ResultStatus>(response.status());
solve_log_ = solution->solve_log;
result_status_ = static_cast<MPSolver::ResultStatus>(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<MPSolutionResponse> 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<MPSolutionResponse> PdlpInterface::DirectlySolveProto(
}
const absl::optional<LazyMutableCopy<MPModelProto>> 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<MPSolverResponseStatus>(response.status()) << " ("
<< response.status() << "): " << response.status_str();
<< ProtoEnumToString<MPSolverResponseStatus>(error_response.status())
<< " (" << error_response.status()
<< "): " << error_response.status_str();
return absl::nullopt;
}
absl::StatusOr<pdlp::SolutionResponseAndLog> solution =
absl::StatusOr<MPSolutionResponse> 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(); }

16
ortools/pdlp/BUILD.bazel
View File

@@ -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"],

7
ortools/pdlp/primal_dual_hybrid_gradient.cc
View File

@@ -1886,7 +1886,7 @@ SolverResult PrimalDualHybridGradient(
std::move(iteration_stats_callback));
}
absl::StatusOr<SolutionResponseAndLog> SolveMpModelProto(
absl::StatusOr<MPSolutionResponse> SolveMpModelProto(
const MPModelProto& model, const PrimalDualHybridGradientParams& params,
const bool relax_integer_variables,
IterationStatsCallback iteration_stats_callback) {
@@ -1944,8 +1944,9 @@ absl::StatusOr<SolutionResponseAndLog> 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 {

20
ortools/pdlp/primal_dual_hybrid_gradient.h
View File

@@ -129,22 +129,16 @@ SolverResult PrimalDualHybridGradient(
std::function<void(const IterationCallbackInfo&)> 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<SolutionResponseAndLog> 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<MPSolutionResponse> SolveMpModelProto(
const MPModelProto& model, const PrimalDualHybridGradientParams& params,
bool relax_integer_variables = false,
std::function<void(const IterationCallbackInfo&)> iteration_stats_callback =

96
ortools/pdlp/primal_dual_hybrid_gradient_test.cc
View File

@@ -1126,24 +1126,26 @@ TEST(SolveProtoTest, SolvesEasyMinimizationLp) {
)pb");
// Using default convergence tolerances and 1 thread.
PrimalDualHybridGradientParams params;
absl::StatusOr<SolutionResponseAndLog> result =
absl::StatusOr<MPSolutionResponse> 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<SolutionResponseAndLog> result =
absl::StatusOr<MPSolutionResponse> 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<SolutionResponseAndLog> result =
absl::StatusOr<MPSolutionResponse> 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<SolutionResponseAndLog> result =
absl::StatusOr<MPSolutionResponse> 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<SolutionResponseAndLog> result =
absl::StatusOr<MPSolutionResponse> 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);

1
ortools/pdlp/quadratic_program_io.cc
View File

@@ -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"

15
ortools/pdlp/sharded_optimization_utils.cc
View File

@@ -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) {

16
ortools/pdlp/sharded_optimization_utils.h
View File

@@ -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.,

289
ortools/pdlp/trust_region.cc
View File

@@ -21,7 +21,6 @@
#include <vector>
#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<double>::infinity()
: 0.0;
}
double UpperBound(int64_t index) const {
return std::isfinite(qp_.constraint_lower_bounds[index])
? std::numeric_limits<double>::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 <typename TrustRegionProblem>
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 <typename TrustRegionProblem>
double CriticalStepSize(const TrustRegionProblem& problem,
const int64_t index) {
if (problem.Objective(index) == 0.0) {
return std::numeric_limits<double>::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 <typename TrustRegionProblem>
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 <typename TrustRegionProblem>
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 <typename ArrayType, typename ValueFunction>
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<std::vector<int64_t>> 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 <typename TrustRegionProblem>
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<int64_t>& 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 <typename TrustRegionProblem>
double RadiusSquaredAtUndecidedComponents(
const TrustRegionProblem& problem,
const std::vector<std::vector<int64_t>>& 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<std::vector<int64_t>>& undecided_components_by_shard) {
return sharder.ParallelSumOverShards([&](const Sharder::Shard& shard) {
const std::vector<int64_t>& 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 <typename TrustRegionProblem>
double RemoveCriticalStepsAbovePivot(
const double pivot, const TrustRegionProblem& problem,
double RemoveCriticalStepsAboveThreshold(
const TrustRegionProblem& problem, const double step_size_threshold,
const Sharder& sharder,
std::vector<std::vector<int64_t>>& undecided_components_by_shard) {
return sharder.ParallelSumOverShards([&](const Sharder::Shard& shard) {
std::vector<int64_t>& 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 <typename TrustRegionProblem>
double RemoveCriticalStepsBelowPivot(
const double pivot, const TrustRegionProblem& problem,
double RemoveCriticalStepsBelowThreshold(
const TrustRegionProblem& problem, const double step_size_threshold,
const Sharder& sharder,
std::vector<std::vector<int64_t>>& undecided_components_by_shard) {
return sharder.ParallelSumOverShards([&](const Sharder::Shard& shard) {
std::vector<int64_t>& 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,

251
ortools/pdlp/trust_region.h
View File

@@ -14,7 +14,16 @@
#ifndef PDLP_TRUST_REGION_H_
#define PDLP_TRUST_REGION_H_
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <limits>
#include <vector>
#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 <typename TrustRegionProblem>
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 <typename TrustRegionProblem>
double CriticalStepSize(const TrustRegionProblem& problem,
const int64_t index) {
if (problem.Objective(index) == 0.0) {
return std::numeric_limits<double>::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 <typename TrustRegionProblem>
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 <typename ArrayType, typename ValueFunction>
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 <typename TrustRegionProblem>
double ComputeInitialUndecidedComponents(
const TrustRegionProblem& problem, int64_t start_index, int64_t end_index,
std::vector<int64_t>& 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 <typename TrustRegionProblem>
double RadiusSquaredOfUndecidedComponents(
const TrustRegionProblem& problem, const double step_size,
const std::vector<int64_t>& 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 <typename TrustRegionProblem>
double RemoveCriticalStepsAboveThreshold(
const TrustRegionProblem& problem, const double step_size_threshold,
std::vector<int64_t>& 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 <typename TrustRegionProblem>
double RemoveCriticalStepsBelowThreshold(
const TrustRegionProblem& problem, const double step_size_threshold,
std::vector<int64_t>& 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<double>::infinity()
: 0.0;
}
double UpperBound(int64_t index) const {
return std::isfinite(qp_.constraint_lower_bounds[index])
? std::numeric_limits<double>::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_

68
ortools/pdlp/trust_region_test.cc
View File

@@ -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;
}
}

196
ortools/sat/encoding.cc
View File

@@ -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<EncodingNode*>& nodes,
SatSolver* solver,
std::deque<EncodingNode>* repository) {
EncodingNode* LazyMergeAllNodeWithPQAndIncreaseLb(
Coefficient weight, const std::vector<EncodingNode*>& nodes,
SatSolver* solver, std::deque<EncodingNode>* repository) {
std::priority_queue<EncodingNode*, std::vector<EncodingNode*>,
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<EncodingNode*>& 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<EncodingNode*> 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<Literal> ReduceNodesAndExtractAssumptions(
@@ -402,14 +453,7 @@ std::vector<Literal> 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<Literal> 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<Literal> ReduceNodesAndExtractAssumptions(
std::vector<Literal> 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<EncodingNode*>& 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<Literal>& core, Coefficient min_weight,
std::vector<EncodingNode*>* 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<Literal>& 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<Literal>& 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<Literal>& core,
Coefficient min_weight,
std::deque<EncodingNode>* repository,
std::vector<EncodingNode*>* 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<EncodingNode*> new_nodes;
std::vector<EncodingNode*> 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

89
ortools/sat/encoding.h
View File

@@ -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<Literal> 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<T> is 32
static_assert(sizeof(EncodingNode) == 72,
"ERROR_EncodingNode_is_not_well_compacted");
#elif defined(_GLIBCXX_DEBUG)
// In debug GNU libstdc++ std::Vector<T> 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<EncodingNode>* repository);
// Same as MergeAllNodesWithDeque() but use a priority queue to merge in
// priority nodes with smaller sizes.
EncodingNode* LazyMergeAllNodeWithPQ(const std::vector<EncodingNode*>& nodes,
SatSolver* solver,
std::deque<EncodingNode>* 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<EncodingNode*>& nodes,
SatSolver* solver, std::deque<EncodingNode>* 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<Literal>& core, Coefficient min_weight,
std::deque<EncodingNode>* repository,
std::vector<EncodingNode*>* 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<Literal>& core,
Coefficient min_weight,
std::deque<EncodingNode>* repository,
std::vector<EncodingNode*>* nodes,
SatSolver* solver);
} // namespace sat
} // namespace operations_research

3
ortools/sat/optimization.cc
View File

@@ -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;

96
ortools/sat/python/swig_helper.cc Normal file
View File

@@ -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 <string>
#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_<SolutionCallback, PySolutionCallback>(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_<SolveWrapper>(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_<CpSatHelper>(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"));
}

6
ortools/sat/sat_solver.cc
View File

@@ -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();

45
ortools/util/python/sorted_interval_list.cc Normal file
View File

@@ -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 <cstdint>
#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_<Domain>(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<int64_t, int64_t>())
.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);
}