diff --git a/CMakeLists.txt b/CMakeLists.txt index 97da2f5b3a..36214948e4 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -336,7 +336,7 @@ include(dotnet) # Since samples mix all languages we must parse them once we have included all # .cmake files -foreach(SAMPLES IN ITEMS algorithms graph glop constraint_solver linear_solver sat) +foreach(SAMPLES IN ITEMS algorithms graph glop constraint_solver linear_solver pdlp sat) add_subdirectory(ortools/${SAMPLES}/samples) endforeach() diff --git a/ortools/pdlp/samples/BUILD.bazel b/ortools/pdlp/samples/BUILD.bazel new file mode 100644 index 0000000000..99968ee5ac --- /dev/null +++ b/ortools/pdlp/samples/BUILD.bazel @@ -0,0 +1,16 @@ +# Copyright 2010-2022 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(":code_samples.bzl", "code_sample_cc") + +code_sample_cc(name = "simple_pdlp_program") diff --git a/ortools/pdlp/samples/CMakeLists.txt b/ortools/pdlp/samples/CMakeLists.txt new file mode 100644 index 0000000000..a3003cb66d --- /dev/null +++ b/ortools/pdlp/samples/CMakeLists.txt @@ -0,0 +1,23 @@ +# Copyright 2010-2022 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. + +if(NOT BUILD_SAMPLES) + return() +endif() + +if(BUILD_CXX_SAMPLES) + file(GLOB CXX_SRCS "*.cc") + foreach(SAMPLE IN LISTS CXX_SRCS) + add_cxx_sample(${SAMPLE}) + endforeach() +endif() diff --git a/ortools/pdlp/samples/code_samples.bzl b/ortools/pdlp/samples/code_samples.bzl new file mode 100644 index 0000000000..17626f6e74 --- /dev/null +++ b/ortools/pdlp/samples/code_samples.bzl @@ -0,0 +1,45 @@ +# Copyright 2010-2022 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. + +"""Helper macro to compile and test code samples.""" + +def code_sample_cc(name): + native.cc_binary( + name = name + "_cc", + srcs = [name + ".cc"], + deps = [ + "//ortools/base", + "//ortools/pdlp:iteration_stats", + "//ortools/pdlp:primal_dual_hybrid_gradient", + "//ortools/pdlp:quadratic_program", + "//ortools/pdlp:solve_log_cc_proto", + "//ortools/pdlp:solvers_cc_proto", + "@eigen//:eigen3", + ], + ) + + native.cc_test( + name = name + "_cc_test", + size = "small", + srcs = [name + ".cc"], + deps = [ + ":" + name + "_cc", + "//ortools/base", + "//ortools/pdlp:iteration_stats", + "//ortools/pdlp:primal_dual_hybrid_gradient", + "//ortools/pdlp:quadratic_program", + "//ortools/pdlp:solve_log_cc_proto", + "//ortools/pdlp:solvers_cc_proto", + "@eigen//:eigen3", + ], + ) diff --git a/ortools/pdlp/samples/simple_pdlp_program.cc b/ortools/pdlp/samples/simple_pdlp_program.cc new file mode 100644 index 0000000000..4a7dc764e7 --- /dev/null +++ b/ortools/pdlp/samples/simple_pdlp_program.cc @@ -0,0 +1,119 @@ +// Copyright 2010-2022 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. + +// Solves a simple LP using PDLP's direct C++ API. +// +// Note: The direct API is generally for advanced use cases. It is matrix-based, +// that is, you specify the LP using matrices and vectors instead of algebraic +// expressions. You can also use PDLP via the algebraic MPSolver API (see +// linear_solver/samples/simple_lp_program.cc). +#include +#include +#include +#include +#include + +#include "Eigen/Core" +#include "Eigen/SparseCore" +#include "ortools/base/init_google.h" +#include "ortools/pdlp/iteration_stats.h" +#include "ortools/pdlp/primal_dual_hybrid_gradient.h" +#include "ortools/pdlp/quadratic_program.h" +#include "ortools/pdlp/solve_log.pb.h" +#include "ortools/pdlp/solvers.pb.h" + +namespace pdlp = ::operations_research::pdlp; + +constexpr double kInfinity = std::numeric_limits::infinity(); + +// Returns a small LP: +// min 5.5 x_0 - 2 x_1 - x_2 + x_3 - 14 s.t. +// 2 x_0 + x_1 + x_2 + 2 x_3 = 12 +// x_0 + x_2 <= 7 +// 4 x_0 >= -4 +// -1 <= 1.5 x_2 - x_3 <= 1 +// -infinity <= x_0 <= infinity +// -2 <= x_1 <= infinity +// -infinity <= x_2 <= 6 +// 2.5 <= x_3 <= 3.5 +pdlp::QuadraticProgram SimpleLp() { + pdlp::QuadraticProgram lp(4, 4); + // "<<" is Eigen's syntax for initialization. + lp.constraint_lower_bounds << 12, -kInfinity, -4, -1; + lp.constraint_upper_bounds << 12, 7, kInfinity, 1; + lp.variable_lower_bounds << -kInfinity, -2, -kInfinity, 2.5; + lp.variable_upper_bounds << kInfinity, kInfinity, 6, 3.5; + const std::vector> + constraint_matrix_triplets = {{0, 0, 2}, {0, 1, 1}, {0, 2, 1}, + {0, 3, 2}, {1, 0, 1}, {1, 2, 1}, + {2, 0, 4}, {3, 2, 1.5}, {3, 3, -1}}; + lp.constraint_matrix.setFromTriplets(constraint_matrix_triplets.begin(), + constraint_matrix_triplets.end()); + lp.objective_vector << 5.5, -2, -1, 1; + lp.objective_offset = -14; + return lp; +} + +int main(int argc, char* argv[]) { + InitGoogle(argv[0], &argc, &argv, /*remove_flags=*/true); + + pdlp::PrimalDualHybridGradientParams params; + // Below are some common parameters to modify. Here, we just re-assign the + // defaults. + params.mutable_termination_criteria() + ->mutable_simple_optimality_criteria() + ->set_eps_optimal_relative(1.0e-6); + params.mutable_termination_criteria() + ->mutable_simple_optimality_criteria() + ->set_eps_optimal_absolute(1.0e-6); + params.mutable_termination_criteria()->set_time_sec_limit(kInfinity); + params.set_num_threads(1); + params.set_verbosity_level(0); + params.mutable_presolve_options()->set_use_glop(false); + + const pdlp::SolverResult result = + pdlp::PrimalDualHybridGradient(SimpleLp(), params); + const pdlp::SolveLog& solve_log = result.solve_log; + + if (solve_log.termination_reason() == pdlp::TERMINATION_REASON_OPTIMAL) { + std::cout << "Solve successful" << std::endl; + } else { + std::cout << "Solve not successful. Status: " + << pdlp::TerminationReason_Name(solve_log.termination_reason()) + << std::endl; + } + + // Solutions vectors are always returned. *However*, their interpretation + // depends on termination_reason! See primal_dual_hybrid_gradient.h for more + // details on what the vectors mean if termination_reason is not + // TERMINATION_REASON_OPTIMAL. + std::cout << "Primal solution:\n" << result.primal_solution << std::endl; + std::cout << "Dual solution:\n" << result.dual_solution << std::endl; + std::cout << "Reduced costs:\n" << result.reduced_costs << std::endl; + + const pdlp::PointType solution_type = solve_log.solution_type(); + std::cout << "Solution type: " << pdlp::PointType_Name(solution_type) + << std::endl; + const std::optional ci = + pdlp::GetConvergenceInformation(solve_log.solution_stats(), + solution_type); + if (ci.has_value()) { + std::cout << "Primal objective: " << ci->primal_objective() << std::endl; + std::cout << "Dual objective: " << ci->dual_objective() << std::endl; + } + + std::cout << "Iterations: " << solve_log.iteration_count() << std::endl; + std::cout << "Solve time (sec): " << solve_log.solve_time_sec() << std::endl; + + return 0; +} diff --git a/ortools/pdlp/samples/simple_pdlp_program.py b/ortools/pdlp/samples/simple_pdlp_program.py new file mode 100644 index 0000000000..74274f70fb --- /dev/null +++ b/ortools/pdlp/samples/simple_pdlp_program.py @@ -0,0 +1,110 @@ +#!/usr/bin/env python3 +# Copyright 2010-2022 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. +"""Solves a simple LP using PDLP's direct Python API. + +Note: The direct API is generally for advanced use cases. It is matrix-based, +that is, you specify the LP using matrices and vectors instead of algebraic +expressions. You can also use PDLP via the algebraic pywraplp API (see +linear_solver/samples/simple_lp_program.py). +""" + +import numpy as np +import scipy.sparse + +from ortools.pdlp import solve_log_pb2 +from ortools.pdlp import solvers_pb2 +from ortools.pdlp.python import pywrap_pdlp +from ortools.init import pywrapinit + + +def simple_lp() -> pywrap_pdlp.QuadraticProgram: + """Returns a small LP. + + min 5.5 x_0 - 2 x_1 - x_2 + x_3 - 14 s.t. + 2 x_0 + x_1 + x_2 + 2 x_3 = 12 + x_0 + x_2 <= 7 + 4 x_0 >= -4 + -1 <= 1.5 x_2 - x_3 <= 1 + -infinity <= x_0 <= infinity + -2 <= x_1 <= infinity + -infinity <= x_2 <= 6 + 2.5 <= x_3 <= 3.5 + """ + lp = pywrap_pdlp.QuadraticProgram() + lp.objective_offset = -14 + lp.objective_vector = [5.5, -2, -1, 1] + lp.constraint_lower_bounds = [12, -np.inf, -4, -1] + lp.constraint_upper_bounds = [12, 7, np.inf, 1] + lp.variable_lower_bounds = [-np.inf, -2, -np.inf, 2.5] + lp.variable_upper_bounds = [np.inf, np.inf, 6, 3.5] + # Most use cases should initialize the sparse constraint matrix without + # constructing a dense matrix first! We use a np.array here for convenience + # only. + constraint_matrix = np.array([[2, 1, 1, 2], [1, 0, 1, 0], [4, 0, 0, 0], + [0, 0, 1.5, -1]]) + lp.constraint_matrix = scipy.sparse.csc_matrix(constraint_matrix) + return lp + + +def main() -> None: + params = solvers_pb2.PrimalDualHybridGradientParams() + # Below are some common parameters to modify. Here, we just re-assign the + # defaults. + optimality_criteria = params.termination_criteria.simple_optimality_criteria + optimality_criteria.eps_optimal_relative = 1.0e-6 + optimality_criteria.eps_optimal_absolute = 1.0e-6 + params.termination_criteria.time_sec_limit = np.inf + params.num_threads = 1 + params.verbosity_level = 0 + params.presolve_options.use_glop = False + + # Call the main solve function. Note that a quirk of the pywrap11 API forces + # us to serialize the `params` and deserialize the `solve_log` proto messages. + result = pywrap_pdlp.primal_dual_hybrid_gradient(simple_lp(), + params.SerializeToString()) + solve_log = solve_log_pb2.SolveLog.FromString(result.solve_log_str) + + if solve_log.termination_reason == solve_log_pb2.TERMINATION_REASON_OPTIMAL: + print('Solve successful') + else: + print( + 'Solve not successful. Status:', + solve_log_pb2.TerminationReason.Name(solve_log.termination_reason)) + + # Solutions vectors are always returned. *However*, their interpretation + # depends on termination_reason! See primal_dual_hybrid_gradient.h for more + # details on what the vectors mean if termination_reason is not + # TERMINATION_REASON_OPTIMAL. + print('Primal solution:', result.primal_solution) + print('Dual solution:', result.dual_solution) + print('Reduced costs:', result.reduced_costs) + + solution_type = solve_log.solution_type + print('Solution type:', solve_log_pb2.PointType.Name(solution_type)) + for ci in solve_log.solution_stats.convergence_information: + if ci.candidate_type == solution_type: + print('Primal objective:', ci.primal_objective) + print('Dual objective:', ci.dual_objective) + + print('Iterations:', solve_log.iteration_count) + print('Solve time (sec):', solve_log.solve_time_sec) + + +if __name__ == '__main__': + pywrapinit.CppBridge.InitLogging('simple_pdlp_program.py') + cpp_flags = pywrapinit.CppFlags() + cpp_flags.logtostderr = True + cpp_flags.log_prefix = False + pywrapinit.CppBridge.SetFlags(cpp_flags) + main()