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

pdlp: add samples

Browse Source
This commit is contained in:
Corentin Le Molgat
2022-10-07 18:24:17 +02:00
parent cecd0b994e
commit 1d5805f8c7
6 changed files with 314 additions and 1 deletions
Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@@ -336,7 +336,7 @@ include(dotnet)
# Since samples mix all languages we must parse them once we have included all
# <language>.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()

16
ortools/pdlp/samples/BUILD.bazel Normal file
View File

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

23
ortools/pdlp/samples/CMakeLists.txt Normal file
View File

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

45
ortools/pdlp/samples/code_samples.bzl Normal file
View File

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

119
ortools/pdlp/samples/simple_pdlp_program.cc Normal file
View File

@@ -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 <cstdint>
#include <iostream>
#include <limits>
#include <optional>
#include <vector>
#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<double>::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<Eigen::Triplet<double, int64_t>>
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<pdlp::ConvergenceInformation> 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;
}

110
ortools/pdlp/samples/simple_pdlp_program.py Normal file
View File

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