From 1d0ff09af8fcc719bb160f678dedcb6dc3f4f716 Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Wed, 9 Jul 2025 14:23:53 +0200 Subject: [PATCH] ortools: backport from main branch --- ortools/algorithms/BUILD.bazel | 1 + ortools/algorithms/hungarian.cc | 1 - ortools/algorithms/hungarian.h | 16 --- ortools/algorithms/knapsack_solver_test.cc | 9 +- ortools/base/proto_enum_utils.h | 1 + ortools/bop/bop_portfolio.h | 1 - ortools/glop/preprocessor.h | 1 - ortools/glop/revised_simplex.cc | 1 - ortools/graph/README.md | 109 +++++++++----------- ortools/graph/cliques.h | 1 - ortools/graph/linear_assignment.h | 1 - ortools/graph/solve_flow_model.cc | 1 - ortools/lp_data/BUILD.bazel | 11 -- ortools/lp_data/README.md | 86 +++++++++++++++ ortools/lp_data/lp_data.h | 1 - ortools/lp_data/sparse.h | 1 - ortools/pdlp/BUILD.bazel | 1 + ortools/pdlp/primal_dual_hybrid_gradient.cc | 97 ++++++++++++----- ortools/python/README.md | 2 +- ortools/util/bitset.h | 3 + ortools/util/permutation.h | 1 - ortools/util/sorted_interval_list.h | 4 +- 22 files changed, 218 insertions(+), 132 deletions(-) diff --git a/ortools/algorithms/BUILD.bazel b/ortools/algorithms/BUILD.bazel index 4f4def32e7..ed81524193 100644 --- a/ortools/algorithms/BUILD.bazel +++ b/ortools/algorithms/BUILD.bazel @@ -240,6 +240,7 @@ cc_test( "//ortools/base:gmock_main", "//ortools/util:time_limit", "@abseil-cpp//absl/base:core_headers", + "@abseil-cpp//absl/types:span", ], ) diff --git a/ortools/algorithms/hungarian.cc b/ortools/algorithms/hungarian.cc index f39d420cc3..b0b2fa109a 100644 --- a/ortools/algorithms/hungarian.cc +++ b/ortools/algorithms/hungarian.cc @@ -16,7 +16,6 @@ #include #include #include -#include #include #include diff --git a/ortools/algorithms/hungarian.h b/ortools/algorithms/hungarian.h index 538027b0c9..5ae49fbbae 100644 --- a/ortools/algorithms/hungarian.h +++ b/ortools/algorithms/hungarian.h @@ -11,18 +11,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -// -// IMPORTANT NOTE: we advise using the code in -// graph/linear_assignment.h whose complexity is -// usually much smaller. -// TODO(user): base this code on LinearSumAssignment. -// -// For each of the four functions declared in this file, in case the input -// parameter 'cost' contains NaN, the function will return without invoking the -// Hungarian algorithm, and the output parameters 'direct_assignment' and -// 'reverse_assignment' will be left unchanged. -// - // An O(n^4) implementation of the Kuhn-Munkres algorithm (a.k.a. the // Hungarian algorithm) for solving the assignment problem. // The assignment problem takes a set of agents, a set of tasks and a @@ -30,10 +18,6 @@ // an optimal (i.e., least cost) assignment of agents to tasks. // The code also enables computing a maximum assignment by changing the // input matrix. -// -// This code is based on (read: translated from) the Java version -// (read: translated from) the Python version at -// http://www.clapper.org/software/python/munkres/. #ifndef OR_TOOLS_ALGORITHMS_HUNGARIAN_H_ #define OR_TOOLS_ALGORITHMS_HUNGARIAN_H_ diff --git a/ortools/algorithms/knapsack_solver_test.cc b/ortools/algorithms/knapsack_solver_test.cc index 1589f20e43..2c7572aead 100644 --- a/ortools/algorithms/knapsack_solver_test.cc +++ b/ortools/algorithms/knapsack_solver_test.cc @@ -18,6 +18,7 @@ #include #include "absl/base/macros.h" +#include "absl/types/span.h" #include "gtest/gtest.h" #include "ortools/util/time_limit.h" @@ -26,8 +27,8 @@ namespace { const int kInvalidSolution = -1; -bool IsSolutionValid(const std::vector& profits, - const std::vector >& weights, +bool IsSolutionValid(absl::Span profits, + absl::Span> weights, const std::vector& capacities, const std::vector& best_solution, int64_t optimal_profit) { @@ -59,7 +60,7 @@ int64_t SolveKnapsackProblemUsingSpecificSolverAndReduction( std::vector profits(profit_array, profit_array + number_of_items); std::vector capacities(capacity_array, capacity_array + number_of_dimensions); - std::vector > weights; + std::vector> weights; for (int i = 0; i < number_of_dimensions; ++i) { const int64_t* one_dimension = weight_array + number_of_items * i; std::vector weights_one_dimension(one_dimension, @@ -484,7 +485,7 @@ TEST(KnapsackSolverTest, SolveTwoDimensionsSettingPrimaryPropagator) { std::vector profits(kProfitArray, kProfitArray + kArraySize); std::vector capacities(kCapacityArray, kCapacityArray + kNumberOfDimensions); - std::vector > weights; + std::vector> weights; for (int i = 0; i < kNumberOfDimensions; ++i) { const int64_t* one_dimension = kWeightArray + kArraySize * i; std::vector weights_one_dimension(one_dimension, diff --git a/ortools/base/proto_enum_utils.h b/ortools/base/proto_enum_utils.h index bdf0331056..4c7dc30ba8 100644 --- a/ortools/base/proto_enum_utils.h +++ b/ortools/base/proto_enum_utils.h @@ -26,6 +26,7 @@ // } // +#include #include #include "google/protobuf/descriptor.pb.h" diff --git a/ortools/bop/bop_portfolio.h b/ortools/bop/bop_portfolio.h index 16f78a616f..085ddaaf9b 100644 --- a/ortools/bop/bop_portfolio.h +++ b/ortools/bop/bop_portfolio.h @@ -56,7 +56,6 @@ class OptimizerSelector; // - LP_FIRST_SOLUTION // - OBJECTIVE_FIRST_SOLUTION // - USER_GUIDED_FIRST_SOLUTION -// - FEASIBILITY_PUMP_FIRST_SOLUTION // - RANDOM_CONSTRAINT_LNS_GUIDED_BY_LP // - RANDOM_VARIABLE_LNS_GUIDED_BY_LP // - RELATION_GRAPH_LNS diff --git a/ortools/glop/preprocessor.h b/ortools/glop/preprocessor.h index dc171cb91d..4e77c5cdfe 100644 --- a/ortools/glop/preprocessor.h +++ b/ortools/glop/preprocessor.h @@ -11,7 +11,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -// // This file contains the presolving code for a LinearProgram. // // A classical reference is: diff --git a/ortools/glop/revised_simplex.cc b/ortools/glop/revised_simplex.cc index fcfe6ef92e..f1aeae58aa 100644 --- a/ortools/glop/revised_simplex.cc +++ b/ortools/glop/revised_simplex.cc @@ -18,7 +18,6 @@ #include #include #include -#include #include #include #include diff --git a/ortools/graph/README.md b/ortools/graph/README.md index 4831200355..7a070bed18 100644 --- a/ortools/graph/README.md +++ b/ortools/graph/README.md @@ -5,76 +5,64 @@ flow problems. It contains in particular: -* well-tuned algorithms (for example, shortest paths and - [Hamiltonian paths](https://en.wikipedia.org/wiki/Hamiltonian_path)). -* hard-to-find algorithms (Hamiltonian paths, push-relabel flow algorithms). -* other, more common algorithms, that are useful to use with graphs from - `util/graph`. +* well-tuned algorithms (for example, shortest paths and + [Hamiltonian paths](https://en.wikipedia.org/wiki/Hamiltonian_path)). +* hard-to-find algorithms (Hamiltonian paths, push-relabel flow algorithms). +* other, more common algorithms, that are useful to use with graphs from + `util/graph`. Generic algorithms for shortest paths: -* [`bounded_dijkstra.h`][bounded_dijkstra_h]: entry point for shortest paths. - This is the preferred implementation for most needs. - -* [`bidirectional_dijkstra.h`][bidirectional_dijkstra_h]: for large graphs, - this implementation might be faster than `bounded_dijkstra.h`. - -* [`shortest_paths.h`][shortest_paths_h]: shortest paths that are computed in - parallel (only if there are several sources). Includes - [Dijkstra](https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm) and - [Bellman-Ford](https://en.wikipedia.org/wiki/Bellman%E2%80%93Ford_algorithm) - algorithms. *Although its name is very generic, only use this implementation - if parallelism makes sense.* +* [`bounded_dijkstra.h`][bounded_dijkstra_h]: entry point for shortest paths. + This is the preferred implementation for most needs. +* [`bidirectional_dijkstra.h`][bidirectional_dijkstra_h]: for large graphs, + this implementation might be faster than `bounded_dijkstra.h`. +* [`shortest_paths.h`][shortest_paths_h]: shortest paths that are computed in + parallel (only if there are several sources). Includes + [Dijkstra](https://en.wikipedia.org/wiki/Dijkstra%27s_algorithm) and + [Bellman-Ford](https://en.wikipedia.org/wiki/Bellman%E2%80%93Ford_algorithm) + algorithms. *Although its name is very generic, only use this implementation + if parallelism makes sense.* Specific algorithms for paths: -* [`dag_shortest_path.h`][dag_shortest_path_h]: shortest paths on directed - acyclic graphs. If you have such a graph, this implementation is likely to - be the fastest. Unlike most implementations, these algorithms have two - interfaces: a "simple" one (list of edges and weights) and a standard one - (taking as input a graph data structure from - [`//ortools/graph/graph.h`][graph_h]). - -* [`dag_constrained_shortest_path.`][dag_constrained_shortest_path_h]: - shortest paths on directed acyclic graphs with resource constraints. - -* [`hamiltonian_path.h`][hamiltonian_path_h]: entry point for computing - minimum [Hamiltonian paths](https://en.wikipedia.org/wiki/Hamiltonian_path) - and cycles on directed graphs with costs on arcs, using a - dynamic-programming algorithm. - -* [`eulerian_path.h`][eulerian_path_h]: entry point for computing minimum - [Eulerian paths](https://en.wikipedia.org/wiki/Eulerian_path) and cycles on - undirected graphs. +* [`dag_shortest_path.h`][dag_shortest_path_h]: shortest paths on directed + acyclic graphs. If you have such a graph, this implementation is likely to be + the fastest. Unlike most implementations, these algorithms have two interfaces + : a "simple" one (list of edges and weights) and a standard one (taking as + input a graph data structure from [`//ortools/graph/graph.h`][graph_h]). +* [`dag_constrained_shortest_path.`][dag_constrained_shortest_path_h]: shortest + paths on directed acyclic graphs with resource constraints. +* [`hamiltonian_path.h`][hamiltonian_path_h]: entry point for computing minimum + [Hamiltonian paths](https://en.wikipedia.org/wiki/Hamiltonian_path) and cycles + on directed graphs with costs on arcs, using a dynamic-programming algorithm. +* [`eulerian_path.h`][eulerian_path_h]: entry point for computing minimum + [Eulerian paths](https://en.wikipedia.org/wiki/Eulerian_path) and cycles on + undirected graphs. Graph decompositions: * [`connected_components.h`][connected_components_h]: entry point for computing connected components in an undirected graph. (It does not need an explicit graph class.) - * [`strongly_connected_components.h`][strongly_connected_components_h]: entry point for computing the strongly connected components of a directed graph. - -* [`cliques.h`][cliques_h]: entry point for computing maximum cliques and - clique covers in a directed graph, based on the Bron-Kerbosch algorithm.(It - does not need an explicit graph class.) +* [`cliques.h`][cliques_h]: entry point for computing maximum cliques and clique + covers in a directed graph, based on the Bron-Kerbosch algorithm.(It does not + need an explicit graph class.) Flow algorithms: -* [`linear_assignment.h`][linear_assignment_h]: entry point for solving linear - sum assignment problems (classical assignment problems where the total cost - is the sum of the costs of each arc used) on directed graphs with arc costs, - based on the Goldberg-Kennedy push-relabel algorithm. - -* [`max_flow.h`][max_flow_h]: entry point for computing maximum flows on - directed graphs with arc capacities, based on the Goldberg-Tarjan - push-relabel algorithm. - -* [`min_cost_flow.h`][min_cost_flow_h]: entry point for computing minimum-cost - flows on directed graphs with arc capacities, arc costs, and - supplies/demands at nodes, based on the Goldberg-Tarjan push-relabel - algorithm. +* [`linear_assignment.h`][linear_assignment_h]: entry point for solving linear + sum assignment problems (classical assignment problems where the total cost is + the sum of the costs of each arc used) on directed graphs with arc costs, + based on the Goldberg-Kennedy push-relabel algorithm. +* [`max_flow.h`][max_flow_h]: entry point for computing maximum flows on + directed graphs with arc capacities, based on the Goldberg-Tarjan push-relabel + algorithm. +* [`min_cost_flow.h`][min_cost_flow_h]: entry point for computing minimum-cost + flows on directed graphs with arc capacities, arc costs, and supplies/demands + at nodes, based on the Goldberg-Tarjan push-relabel algorithm. ## Class design @@ -201,21 +189,18 @@ node. ## Wrappers -* [`python`](python): the SWIG code that makes the wrapper available in Python - and its unit tests. - -* [`java`](java): the SWIG code that makes the wrapper available in Java and - its unit tests. - -* [`csharp`](csharp): the SWIG code that makes the wrapper available in C# and - its unit tests. +* [`python`](python): This directory contains the `pybind11` bindings to make + these C++ libraries available in Python. +* [`java`](java): the SWIG code that makes the wrapper available in Java and its + unit tests. +* [`csharp`](csharp): the SWIG code that makes the wrapper available in C# and + its unit tests. ## Samples You can find some canonical examples in [`samples`][samples]. - [graph_h]: ../graph/graph.h [bounded_dijkstra_h]: ../graph/bounded_dijkstra.h [bidirectional_dijkstra_h]: ../graph/bidirectional_dijkstra.h diff --git a/ortools/graph/cliques.h b/ortools/graph/cliques.h index 387c49a42c..6afe4a9f8a 100644 --- a/ortools/graph/cliques.h +++ b/ortools/graph/cliques.h @@ -11,7 +11,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -// // Maximal clique algorithms, based on the Bron-Kerbosch algorithm. // See http://en.wikipedia.org/wiki/Bron-Kerbosch_algorithm // and diff --git a/ortools/graph/linear_assignment.h b/ortools/graph/linear_assignment.h index 188a46ddab..c77ad79b1f 100644 --- a/ortools/graph/linear_assignment.h +++ b/ortools/graph/linear_assignment.h @@ -11,7 +11,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -// // An implementation of a cost-scaling push-relabel algorithm for the // assignment problem (minimum-cost perfect bipartite matching), from // the paper of Goldberg and Kennedy (1995). diff --git a/ortools/graph/solve_flow_model.cc b/ortools/graph/solve_flow_model.cc index 31e98c162f..fd42f9fc6e 100644 --- a/ortools/graph/solve_flow_model.cc +++ b/ortools/graph/solve_flow_model.cc @@ -11,7 +11,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -// // This code loads flow-graph models (as Dimacs file or binary FlowModel proto) // and solves them with the OR-tools flow algorithms. // diff --git a/ortools/lp_data/BUILD.bazel b/ortools/lp_data/BUILD.bazel index 6e3ad9bd97..619de54c9a 100644 --- a/ortools/lp_data/BUILD.bazel +++ b/ortools/lp_data/BUILD.bazel @@ -254,17 +254,6 @@ cc_library( ], ) -#cc_library( -# name = "lp_constraint_classifier", -# srcs = ["lp_constraint_classifier.cc"], -# hdrs = ["lp_constraint_classifier.h"], -# copts = SAFE_FP_CODE, -# deps = [ -# ":lp_data", -# "//ortools/util:fp_utils", -# ], -#) - cc_library( name = "lp_print_utils", srcs = ["lp_print_utils.cc"], diff --git a/ortools/lp_data/README.md b/ortools/lp_data/README.md index e69de29bb2..52d4764280 100644 --- a/ortools/lp_data/README.md +++ b/ortools/lp_data/README.md @@ -0,0 +1,86 @@ +# LP Data + +This directory contains a rich collection of C++ libraries for handling Linear +Programming (LP) data structures. + +It provides core components for representing, manipulating, and solving linear +programs, with a focus on efficient handling of sparse data and various utility +functions for pre-processing and analysis. + +## Core Data Structures + +This set of libraries provides the fundamental building blocks for representing +and working with linear programming problems. + +* [`lp_types.h`][lp_types_h]: Defines common types and constants used throughout + the linear programming solver. +* [`lp_data.h`][lp_data_h]: Provides the main `LinearProgram` class for storing + the complete data of a linear program, including the objective function, + constraint matrix, and variable bounds. +* [`lp_utils.h`][lp_utils_h]: Contains basic utility functions for operations on + fractional numbers and row/column vectors. + +## Sparse Data Representation + +Given that large-scale linear programs are often sparse, this directory offers a +suite of libraries for efficient sparse data handling. + +* [`sparse.h`][sparse_h]: Implements data structures for sparse matrices, based + on well-established references in the field of direct methods for sparse + matrices. +* [`sparse_vector.h`][sparse_vector_h]: Provides classes to represent sparse + vectors efficiently. +* [`sparse_column.h`][sparse_column_h] & [`sparse_row.h`][sparse_row_h]: + Specializations of sparse vectors for column-oriented and row-oriented matrix + storage schemes. +* [`scattered_vector.h`][scattered_vector_h]: Implements vectors that offer a + sparse interface to what is internally a dense storage, which can be useful + for certain computations. + +## LP Solvers and Utilities + +A collection of tools for preprocessing, analyzing, and manipulating linear +programs. + +* [`matrix_scaler.h`][matrix_scaler_h]: Provides the `SparseMatrixScaler` class, + which scales a `SparseMatrix` to improve numerical stability during the + solving process. +* [`lp_decomposer.h`][lp_decomposer_h]: Implements a tool to decompose a large + `LinearProgram` into several smaller, independent subproblems by identifying + disconnected components in the constraint matrix. +* [`permutation.h`][permutation_h]: Contains utilities for handling row and + column permutations on LP data structures. + +## Parsers and I/O Utilities + +This group of libraries handles reading and writing LP data in various formats. + +* [`lp_parser.h`][lp_parser_h]: A simple parser for creating a linear program + from a string representation. +* [`mps_reader.h`][mps_reader_h]: A reader for the industry-standard MPS file + format for mathematical programming problems. +* [`sol_reader.h`][sol_reader_h]: A reader for .sol files, which are used to + specify solution values for a given model. +* [`proto_utils.h`][proto_utils_h]: Provides utilities to convert + `LinearProgram` objects to and from the MPModelProto protobuf format. +* [`lp_print_utils.h`][lp_print_utils_h]: Contains utilities to display linear + expressions in a human-readable way, including rational approximations. + + + +[lp_types_h]: ../lp_data/lp_types.h +[lp_data_h]: ../lp_data/lp_data.h +[lp_utils_h]: ../lp_data/lp_utils.h +[sparse_h]: ../lp_data/sparse.h +[sparse_vector_h]: ../lp_data/sparse_vector.h +[sparse_column_h]: ../lp_data/sparse_column.h +[sparse_row_h]: ../lp_data/sparse_row.h +[scattered_vector_h]: ../lp_data/scattered_vector.h +[matrix_scaler_h]: ../lp_data/matrix_scaler.h +[lp_decomposer_h]: ../lp_data/lp_decomposer.h +[permutation_h]: ../lp_data/permutation.h +[lp_parser_h]: ../lp_data/lp_parser.h +[mps_reader_h]: ../lp_data/mps_reader.h +[sol_reader_h]: ../lp_data/sol_reader.h +[proto_utils_h]: ../lp_data/proto_utils.h +[lp_print_utils_h]: ../lp_data/lp_print_utils.h diff --git a/ortools/lp_data/lp_data.h b/ortools/lp_data/lp_data.h index 236423daa1..478a1a35be 100644 --- a/ortools/lp_data/lp_data.h +++ b/ortools/lp_data/lp_data.h @@ -11,7 +11,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -// // Storage classes for Linear Programs. // // LinearProgram stores the complete data for a Linear Program: diff --git a/ortools/lp_data/sparse.h b/ortools/lp_data/sparse.h index fecc42a79c..cc856fc3a8 100644 --- a/ortools/lp_data/sparse.h +++ b/ortools/lp_data/sparse.h @@ -11,7 +11,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -// // The following are very good references for terminology, data structures, // and algorithms: // diff --git a/ortools/pdlp/BUILD.bazel b/ortools/pdlp/BUILD.bazel index 7a19ebc39a..0059d6d122 100644 --- a/ortools/pdlp/BUILD.bazel +++ b/ortools/pdlp/BUILD.bazel @@ -144,6 +144,7 @@ cc_library( "//ortools/lp_data:proto_utils", "//ortools/util:logging", "@abseil-cpp//absl/algorithm:container", + "@abseil-cpp//absl/base:nullability", "@abseil-cpp//absl/status", "@abseil-cpp//absl/status:statusor", "@abseil-cpp//absl/strings", diff --git a/ortools/pdlp/primal_dual_hybrid_gradient.cc b/ortools/pdlp/primal_dual_hybrid_gradient.cc index d166ef7700..b86f3e9f1c 100644 --- a/ortools/pdlp/primal_dual_hybrid_gradient.cc +++ b/ortools/pdlp/primal_dual_hybrid_gradient.cc @@ -53,6 +53,7 @@ #include "Eigen/Core" #include "Eigen/SparseCore" #include "absl/algorithm/container.h" +#include "absl/base/nullability.h" #include "absl/log/check.h" #include "absl/log/log.h" #include "absl/status/status.h" @@ -619,7 +620,8 @@ class Solver { NextSolutionAndDelta ComputeNextDualSolution( double dual_step_size, double extrapolation_factor, - const NextSolutionAndDelta& next_primal) const; + const NextSolutionAndDelta& next_primal_solution, + const VectorXd* absl_nullable next_primal_product = nullptr) const; std::pair ComputeMovementTerms( const VectorXd& delta_primal, const VectorXd& delta_dual) const; @@ -630,6 +632,10 @@ class Solver { double ComputeNonlinearity(const VectorXd& delta_primal, const VectorXd& next_dual_product) const; + // Sets current_primal_product_ and current_dual_product_ based on + // current_primal_solution_ and current_dual_solution_ respectively. + void SetCurrentPrimalAndDualProducts(); + // Creates all the simple-to-compute statistics in stats. IterationStats CreateSimpleIterationStats(RestartChoice restart_used) const; @@ -734,6 +740,9 @@ class Solver { WallTimer timer_; int iterations_completed_; int num_rejected_steps_; + // A cache of `constraint_matrix * current_primal_solution_`. + // Malitsky-Pock linesearch only. + std::optional current_primal_product_; // A cache of `constraint_matrix.transpose() * current_dual_solution_`. VectorXd current_dual_product_; // The primal point at which the algorithm was last restarted from, or @@ -1870,31 +1879,41 @@ Solver::NextSolutionAndDelta Solver::ComputeNextPrimalSolution( Solver::NextSolutionAndDelta Solver::ComputeNextDualSolution( double dual_step_size, double extrapolation_factor, - const NextSolutionAndDelta& next_primal_solution) const { + const NextSolutionAndDelta& next_primal_solution, + const VectorXd* absl_nullable next_primal_product) const { const int64_t dual_size = ShardedWorkingQp().DualSize(); NextSolutionAndDelta result = { .value = VectorXd(dual_size), .delta = VectorXd(dual_size), }; const QuadraticProgram& qp = WorkingQp(); - VectorXd extrapolated_primal(ShardedWorkingQp().PrimalSize()); - ShardedWorkingQp().PrimalSharder().ParallelForEachShard( - [&](const Sharder::Shard& shard) { - shard(extrapolated_primal) = - (shard(next_primal_solution.value) + - extrapolation_factor * shard(next_primal_solution.delta)); - }); - // TODO(user): Refactor this multiplication so that we only do one matrix - // vector multiply for the primal variable. This only applies to Malitsky and - // Pock and not to the adaptive step size rule. + std::optional extrapolated_primal; + if (!next_primal_product) { + extrapolated_primal.emplace(ShardedWorkingQp().PrimalSize()); + ShardedWorkingQp().PrimalSharder().ParallelForEachShard( + [&](const Sharder::Shard& shard) { + shard(*extrapolated_primal) = + (shard(next_primal_solution.value) + + extrapolation_factor * shard(next_primal_solution.delta)); + }); + } ShardedWorkingQp().TransposedConstraintMatrixSharder().ParallelForEachShard( [&](const Sharder::Shard& shard) { - VectorXd temp = - shard(current_dual_solution_) - - dual_step_size * - shard(ShardedWorkingQp().TransposedConstraintMatrix()) - .transpose() * - extrapolated_primal; + VectorXd temp; + if (next_primal_product) { + CHECK(current_primal_product_.has_value()); + temp = shard(current_dual_solution_) - + dual_step_size * + (-extrapolation_factor * shard(*current_primal_product_) + + (extrapolation_factor + 1) * shard(*next_primal_product)); + } else { + temp = shard(current_dual_solution_) - + dual_step_size * + shard(ShardedWorkingQp().TransposedConstraintMatrix()) + .transpose() * + extrapolated_primal.value(); + } + // Each element of the argument of `.cwiseMin()` is the critical point // of the respective 1D minimization problem if it's negative. // Likewise the argument to the `.cwiseMax()` is the critical point if @@ -1937,6 +1956,21 @@ double Solver::ComputeNonlinearity(const VectorXd& delta_primal, }); } +void Solver::SetCurrentPrimalAndDualProducts() { + if (params_.linesearch_rule() == + PrimalDualHybridGradientParams::MALITSKY_POCK_LINESEARCH_RULE) { + current_primal_product_ = TransposedMatrixVectorProduct( + ShardedWorkingQp().TransposedConstraintMatrix(), + current_primal_solution_, + ShardedWorkingQp().TransposedConstraintMatrixSharder()); + } else { + current_primal_product_.reset(); + } + current_dual_product_ = TransposedMatrixVectorProduct( + WorkingQp().constraint_matrix, current_dual_solution_, + ShardedWorkingQp().ConstraintMatrixSharder()); +} + IterationStats Solver::CreateSimpleIterationStats( RestartChoice restart_used) const { IterationStats stats; @@ -1977,8 +2011,10 @@ LocalizedLagrangianBounds Solver::ComputeLocalizedBoundsAtCurrent() const { ShardedWorkingQp(), current_primal_solution_, current_dual_solution_, PrimalDualNorm::kEuclideanNorm, primal_weight_, distance_traveled_by_current, - /*primal_product=*/nullptr, ¤t_dual_product_, - params_.use_diagonal_qp_trust_region_solver(), + /*primal_product=*/current_primal_product_.has_value() + ? ¤t_primal_product_.value() + : nullptr, + ¤t_dual_product_, params_.use_diagonal_qp_trust_region_solver(), params_.diagonal_qp_trust_region_solver_tolerance()); } @@ -2225,9 +2261,7 @@ void Solver::ApplyRestartChoice(const RestartChoice restart_to_apply) { } current_primal_solution_ = primal_average_.ComputeAverage(); current_dual_solution_ = dual_average_.ComputeAverage(); - current_dual_product_ = TransposedMatrixVectorProduct( - WorkingQp().constraint_matrix, current_dual_solution_, - ShardedWorkingQp().ConstraintMatrixSharder()); + SetCurrentPrimalAndDualProducts(); break; } primal_weight_ = ComputeNewPrimalWeight(); @@ -2443,6 +2477,14 @@ InnerStepOutcome Solver::TakeMalitskyPockStep() { params_.malitsky_pock_parameters().linesearch_contraction_factor(); const double dual_weight = primal_weight_ * primal_weight_; int inner_iterations = 0; + VectorXd next_primal_product(current_dual_solution_.size()); + ShardedWorkingQp().TransposedConstraintMatrixSharder().ParallelForEachShard( + [&](const Sharder::Shard& shard) { + shard(next_primal_product) = + shard(ShardedWorkingQp().TransposedConstraintMatrix()).transpose() * + next_primal_solution.value; + }); + for (bool accepted_step = false; !accepted_step; ++inner_iterations) { if (inner_iterations >= 60) { LogInnerIterationLimitHit(); @@ -2454,7 +2496,7 @@ InnerStepOutcome Solver::TakeMalitskyPockStep() { new_primal_step_size / primal_step_size; NextSolutionAndDelta next_dual_solution = ComputeNextDualSolution( dual_weight * new_primal_step_size, new_last_two_step_sizes_ratio, - next_primal_solution); + next_primal_solution, &next_primal_product); VectorXd next_dual_product = TransposedMatrixVectorProduct( WorkingQp().constraint_matrix, next_dual_solution.value, @@ -2482,6 +2524,7 @@ InnerStepOutcome Solver::TakeMalitskyPockStep() { current_primal_solution_ = std::move(next_primal_solution.value); current_dual_solution_ = std::move(next_dual_solution.value); current_dual_product_ = std::move(next_dual_product); + current_primal_product_ = std::move(next_primal_product); primal_average_.Add(current_primal_solution_, /*weight=*/new_primal_step_size); dual_average_.Add(current_dual_solution_, @@ -2555,6 +2598,7 @@ InnerStepOutcome Solver::TakeAdaptiveStep() { current_primal_solution_ = std::move(next_primal_solution.value); current_dual_solution_ = std::move(next_dual_solution.value); current_dual_product_ = std::move(next_dual_product); + current_primal_product_.reset(); current_primal_delta_ = std::move(next_primal_solution.delta); current_dual_delta_ = std::move(next_dual_solution.delta); primal_average_.Add(current_primal_solution_, /*weight=*/step_size_); @@ -2620,6 +2664,7 @@ InnerStepOutcome Solver::TakeConstantSizeStep() { current_primal_solution_ = std::move(next_primal_solution.value); current_dual_solution_ = std::move(next_dual_solution.value); current_dual_product_ = std::move(next_dual_product); + current_primal_product_.reset(); current_primal_delta_ = std::move(next_primal_solution.delta); current_dual_delta_ = std::move(next_dual_solution.delta); primal_average_.Add(current_primal_solution_, /*weight=*/step_size_); @@ -2980,9 +3025,7 @@ SolverResult Solver::Solve(const IterationType iteration_type, // restart. ratio_last_two_step_sizes_ = 1; - current_dual_product_ = TransposedMatrixVectorProduct( - WorkingQp().constraint_matrix, current_dual_solution_, - ShardedWorkingQp().ConstraintMatrixSharder()); + SetCurrentPrimalAndDualProducts(); // This is set to true if we can't proceed any more because of numerical // issues. We may or may not have found the optimal solution. diff --git a/ortools/python/README.md b/ortools/python/README.md index e95202f480..82c903afe9 100644 --- a/ortools/python/README.md +++ b/ortools/python/README.md @@ -18,7 +18,7 @@ This project aim to explain how you build a Python native wheel package using ## Requirement -You'll need "Python >= 3.6" and few python modules ("wheel" and "absl-py"). +You'll need "Python >= 3.9" and few python modules ("wheel" and "absl-py"). ## Directory Layout diff --git a/ortools/util/bitset.h b/ortools/util/bitset.h index 35f3f47e2f..909e0a5551 100644 --- a/ortools/util/bitset.h +++ b/ortools/util/bitset.h @@ -884,6 +884,9 @@ class SparseBitset { // A bit hacky for really hot loop. typename Bitset64::View BitsetView() { return bitset_.view(); } + typename Bitset64::ConstView BitsetConstView() { + return bitset_.const_view(); + } void SetUnsafe(typename Bitset64::View view, IntegerType index) { view.Set(index); to_clear_.push_back(index); diff --git a/ortools/util/permutation.h b/ortools/util/permutation.h index e44906b324..f8cab1b074 100644 --- a/ortools/util/permutation.h +++ b/ortools/util/permutation.h @@ -11,7 +11,6 @@ // See the License for the specific language governing permissions and // limitations under the License. -// // Classes for permuting indexable, ordered containers of data without // depending on that data to be accessible in any particular way. The // client needs to give us two things: diff --git a/ortools/util/sorted_interval_list.h b/ortools/util/sorted_interval_list.h index f07dca7c71..fb62e30d27 100644 --- a/ortools/util/sorted_interval_list.h +++ b/ortools/util/sorted_interval_list.h @@ -724,7 +724,9 @@ class ClosedInterval::Iterator { // arithmetic. uint64_t current_; }; - +#if __cplusplus >= 202002L +static_assert(std::input_iterator); +#endif // begin()/end() are required for iteration over ClosedInterval in a range for // loop. inline ClosedInterval::Iterator begin(ClosedInterval interval) {