Update pdlp

ortools/pdlp/BUILD.bazel

@@ -267,6 +267,7 @@ cc_test(
         ":gtest_main",
         ":solvers_cc_proto",
         ":solvers_proto_validation",
+        "//ortools/base:protobuf_util",
         "@com_google_absl//absl/status",
     ],
 )

ortools/pdlp/iteration_stats.cc

@@ -16,20 +16,15 @@
 #include <algorithm>
 #include <cmath>
 #include <cstdint>
-#include <limits>
 #include <optional>
 #include <random>
-#include <string>
 #include <utility>
 #include <vector>
 
 #include "Eigen/Core"
 #include "Eigen/SparseCore"
 #include "absl/random/distributions.h"
-#include "absl/strings/str_cat.h"
-#include "absl/strings/str_format.h"
-#include "absl/strings/string_view.h"
-#include "ortools/base/logging.h"
+#include "ortools/base/check.h"
 #include "ortools/base/mathutil.h"
 #include "ortools/pdlp/quadratic_program.h"
 #include "ortools/pdlp/sharded_quadratic_program.h"

ortools/pdlp/iteration_stats.h

@@ -14,9 +14,7 @@
 #ifndef PDLP_ITERATION_STATS_H_
 #define PDLP_ITERATION_STATS_H_
 
-#include <limits>
 #include <optional>
-#include <string>
 #include <vector>
 
 #include "Eigen/Core"
@@ -79,17 +77,17 @@ Eigen::VectorXd ReducedCosts(const ShardedQuadraticProgram& scaled_sharded_qp,
                              bool use_zero_primal_objective = false);
 
 // Finds and returns the ConvergenceInformation with the specified
-// candidate_type, or absl::nullopt if no such candidate exists.
+// candidate_type, or std::nullopt if no such candidate exists.
 std::optional<ConvergenceInformation> GetConvergenceInformation(
     const IterationStats& stats, PointType candidate_type);
 
 // Finds and returns the InfeasibilityInformation with the specified
-// candidate_type, or absl::nullopt if no such candidate exists.
+// candidate_type, or std::nullopt if no such candidate exists.
 std::optional<InfeasibilityInformation> GetInfeasibilityInformation(
     const IterationStats& stats, PointType candidate_type);
 
 // Finds and returns the PointMetadata with the specified
-// point_type, or absl::nullopt if no such point exists.
+// point_type, or std::nullopt if no such point exists.
 std::optional<PointMetadata> GetPointMetadata(const IterationStats& stats,
                                               PointType point_type);
 

ortools/pdlp/iteration_stats_test.cc

@@ -14,7 +14,7 @@
 #include "ortools/pdlp/iteration_stats.h"
 
 #include <cmath>
-#include <limits>
+#include <optional>
 #include <utility>
 
 #include "Eigen/Core"
@@ -218,7 +218,7 @@ TEST(GetConvergenceInformation, GetsCorrectEntry) {
 
   EXPECT_THAT(
       GetConvergenceInformation(test_stats, POINT_TYPE_ITERATE_DIFFERENCE),
-      Eq(absl::nullopt));
+      Eq(std::nullopt));
 }
 
 TEST(GetInfeasibilityInformation, GetsCorrectEntry) {
@@ -246,7 +246,7 @@ TEST(GetInfeasibilityInformation, GetsCorrectEntry) {
 
   EXPECT_THAT(
       GetInfeasibilityInformation(test_stats, POINT_TYPE_ITERATE_DIFFERENCE),
-      Eq(absl::nullopt));
+      Eq(std::nullopt));
 }
 
 TEST(GetPointMetadata, GetsCorrectEntry) {
@@ -273,7 +273,7 @@ TEST(GetPointMetadata, GetsCorrectEntry) {
   EXPECT_EQ(current_info->active_primal_variable_count(), 1);
 
   EXPECT_THAT(GetPointMetadata(test_stats, POINT_TYPE_ITERATE_DIFFERENCE),
-              Eq(absl::nullopt));
+              Eq(std::nullopt));
 }
 
 }  // namespace

ortools/pdlp/primal_dual_hybrid_gradient.cc

@@ -174,9 +174,9 @@ std::string ToString(const IterationStats& iter_stats,
   }
   if (convergence_information.has_value()) {
     const RelativeConvergenceInformation relative_information =
-        ComputeRelativeResiduals(termination_criteria.eps_optimal_absolute(),
-                                 termination_criteria.eps_optimal_relative(),
-                                 bound_norms, *convergence_information);
+        ComputeRelativeResiduals(
+            EffectiveOptimalityCriteria(termination_criteria), bound_norms,
+            *convergence_information);
     return absl::StrCat(iteration_string, " | ",
                         ToString(*convergence_information, relative_information,
                                  termination_criteria.optimality_norm()));
@@ -199,9 +199,9 @@ std::string ToShortString(const IterationStats& iter_stats,
   }
   if (convergence_information.has_value()) {
     const RelativeConvergenceInformation relative_information =
-        ComputeRelativeResiduals(termination_criteria.eps_optimal_absolute(),
-                                 termination_criteria.eps_optimal_relative(),
-                                 bound_norms, *convergence_information);
+        ComputeRelativeResiduals(
+            EffectiveOptimalityCriteria(termination_criteria), bound_norms,
+            *convergence_information);
     return absl::StrCat(
         iteration_string, " | ",
         ToShortString(*convergence_information, relative_information,
@@ -1273,7 +1273,7 @@ std::optional<SolverResult> Solver::MajorIterationAndTerminationCheck(
     *solve_log.add_iteration_stats() = stats;
   }
   ApplyRestartChoice(restart);
-  return absl::nullopt;
+  return std::nullopt;
 }
 
 void Solver::ResetAverageToCurrent() {
@@ -1589,19 +1589,19 @@ std::optional<TerminationReason> Solver::ApplyPresolveIfEnabled(
     std::optional<PrimalAndDualSolution>* const initial_solution) {
   const bool presolve_enabled = params_.presolve_options().use_glop();
   if (!presolve_enabled) {
-    return absl::nullopt;
+    return std::nullopt;
   }
   if (!IsLinearProgram(WorkingQp())) {
     LOG(WARNING)
         << "Skipping presolve, which is only supported for linear programs";
-    return absl::nullopt;
+    return std::nullopt;
   }
   absl::StatusOr<MPModelProto> model = QpToMpModelProto(WorkingQp());
   if (!model.ok()) {
     LOG(WARNING)
         << "Skipping presolve because of error converting to MPModelProto: "
         << model.status();
-    return absl::nullopt;
+    return std::nullopt;
   }
   if (initial_solution->has_value()) {
     LOG(WARNING) << "Ignoring initial solution. Initial solutions "
@@ -1643,7 +1643,7 @@ std::optional<TerminationReason> Solver::ApplyPresolveIfEnabled(
     row_scaling_vec_ = OnesVector(sharded_working_qp_.DualSharder());
     return GlopStatusToTerminationReason(presolve_info_->preprocessor.status());
   }
-  return absl::nullopt;
+  return std::nullopt;
 }
 
 PrimalAndDualSolution Solver::RecoverOriginalSolution(
@@ -1933,7 +1933,7 @@ SolverResult Solver::Solve(
     double inverse_step_size;
     const auto lipschitz_result =
         EstimateMaximumSingularValueOfConstraintMatrix(
-            sharded_working_qp_, absl::nullopt, absl::nullopt,
+            sharded_working_qp_, std::nullopt, std::nullopt,
             /*desired_relative_error=*/0.2, /*failure_probability=*/0.0005,
             random);
     // With high probability, the estimate of the lipschitz term is within
@@ -2032,7 +2032,7 @@ SolverResult PrimalDualHybridGradient(
     QuadraticProgram qp, const PrimalDualHybridGradientParams& params,
     const std::atomic<bool>* interrupt_solve,
     IterationStatsCallback iteration_stats_callback) {
-  return PrimalDualHybridGradient(std::move(qp), params, absl::nullopt,
+  return PrimalDualHybridGradient(std::move(qp), params, std::nullopt,
                                   interrupt_solve,
                                   std::move(iteration_stats_callback));
 }

ortools/pdlp/primal_dual_hybrid_gradient.h

@@ -19,8 +19,6 @@
 #include <optional>
 
 #include "Eigen/Core"
-#include "absl/status/statusor.h"
-#include "ortools/linear_solver/linear_solver.pb.h"
 #include "ortools/lp_data/lp_data.h"
 #include "ortools/pdlp/quadratic_program.h"
 #include "ortools/pdlp/solve_log.pb.h"

ortools/pdlp/primal_dual_hybrid_gradient_test.cc

@@ -18,6 +18,7 @@
 #include <cmath>
 #include <cstdint>
 #include <limits>
+#include <optional>
 #include <string>
 #include <tuple>
 #include <utility>
@@ -30,6 +31,7 @@
 #include "absl/strings/str_cat.h"
 #include "gmock/gmock.h"
 #include "gtest/gtest.h"
+#include "ortools/base/check.h"
 #include "ortools/base/logging.h"
 #include "ortools/glop/parameters.pb.h"
 #include "ortools/linear_solver/linear_solver.pb.h"
@@ -75,14 +77,19 @@ PrimalDualHybridGradientParams CreateSolverParams(
         PrimalDualHybridGradientParams::MALITSKY_POCK_LINESEARCH_RULE);
   }
 
-  params.mutable_termination_criteria()->set_eps_optimal_relative(0.0);
+  params.mutable_termination_criteria()
+      ->mutable_simple_optimality_criteria()
+      ->set_eps_optimal_relative(0.0);
   if (use_iteration_limit) {
     // This effectively forces convergence on the iteration limit only.
     params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
-    params.mutable_termination_criteria()->set_eps_optimal_absolute(0.0);
+    params.mutable_termination_criteria()
+        ->mutable_simple_optimality_criteria()
+        ->set_eps_optimal_absolute(0.0);
   } else {
-    params.mutable_termination_criteria()->set_eps_optimal_absolute(
-        eps_optimal_absolute);
+    params.mutable_termination_criteria()
+        ->mutable_simple_optimality_criteria()
+        ->set_eps_optimal_absolute(eps_optimal_absolute);
   }
   if (use_diagonal_qp_trust_region_solver) {
     params.set_use_diagonal_qp_trust_region_solver(true);
@@ -613,8 +620,12 @@ PrimalDualHybridGradientParams ParamsWithNoLimits() {
   PrimalDualHybridGradientParams params;
   // This disables the termination limits. A termination criteria must be set
   // for the solver to terminate.
-  params.mutable_termination_criteria()->set_eps_optimal_relative(0.0);
-  params.mutable_termination_criteria()->set_eps_optimal_absolute(0.0);
+  params.mutable_termination_criteria()
+      ->mutable_simple_optimality_criteria()
+      ->set_eps_optimal_relative(0.0);
+  params.mutable_termination_criteria()
+      ->mutable_simple_optimality_criteria()
+      ->set_eps_optimal_absolute(0.0);
   params.set_record_iteration_stats(true);
   return params;
 }
@@ -882,24 +893,24 @@ TEST(PrimalDualHybridGradientTest,
   EXPECT_EQ(output.solve_log.iteration_stats().size(), iteration_limit + 1);
   for (const auto& stats : output.solve_log.iteration_stats()) {
     EXPECT_NE(GetConvergenceInformation(stats, POINT_TYPE_CURRENT_ITERATE),
-              absl::nullopt);
+              std::nullopt);
     EXPECT_NE(GetInfeasibilityInformation(stats, POINT_TYPE_CURRENT_ITERATE),
-              absl::nullopt);
+              std::nullopt);
     EXPECT_NE(GetPointMetadata(stats, POINT_TYPE_CURRENT_ITERATE),
-              absl::nullopt);
+              std::nullopt);
     if (stats.iteration_number() > 0) {
       EXPECT_NE(GetConvergenceInformation(stats, POINT_TYPE_AVERAGE_ITERATE),
-                absl::nullopt);
+                std::nullopt);
       EXPECT_NE(GetInfeasibilityInformation(stats, POINT_TYPE_AVERAGE_ITERATE),
-                absl::nullopt);
+                std::nullopt);
       EXPECT_NE(GetPointMetadata(stats, POINT_TYPE_AVERAGE_ITERATE),
-                absl::nullopt);
+                std::nullopt);
 
       EXPECT_NE(
           GetInfeasibilityInformation(stats, POINT_TYPE_ITERATE_DIFFERENCE),
-          absl::nullopt);
+          std::nullopt);
       EXPECT_NE(GetPointMetadata(stats, POINT_TYPE_ITERATE_DIFFERENCE),
-                absl::nullopt);
+                std::nullopt);
     }
   }
 }
@@ -1214,7 +1225,9 @@ PrimalAndDualSolution TinyLpSolution() {
 
 TEST(PrimalDualHybridGradientTest, WarmStartedAtOptimum) {
   PrimalDualHybridGradientParams params = ParamsWithNoLimits();
-  params.mutable_termination_criteria()->set_eps_optimal_absolute(1.0e-10);
+  params.mutable_termination_criteria()
+      ->mutable_simple_optimality_criteria()
+      ->set_eps_optimal_absolute(1.0e-10);
 
   SolverResult output =
       PrimalDualHybridGradient(TinyLp(), params, TinyLpSolution());
@@ -1260,8 +1273,12 @@ TEST(PrimalDualHybridGradientTest, EmptyQp) {
 TEST(PrimalDualHybridGradientTest, RespectsInterrupt) {
   std::atomic<bool> interrupt_solve;
   PrimalDualHybridGradientParams params;
-  params.mutable_termination_criteria()->set_eps_optimal_absolute(0.0);
-  params.mutable_termination_criteria()->set_eps_optimal_relative(0.0);
+  params.mutable_termination_criteria()
+      ->mutable_simple_optimality_criteria()
+      ->set_eps_optimal_absolute(0.0);
+  params.mutable_termination_criteria()
+      ->mutable_simple_optimality_criteria()
+      ->set_eps_optimal_relative(0.0);
 
   interrupt_solve.store(true);
   const SolverResult output =
@@ -1273,8 +1290,12 @@ TEST(PrimalDualHybridGradientTest, RespectsInterrupt) {
 TEST(PrimalDualHybridGradientTest, RespectsInterruptFromCallback) {
   std::atomic<bool> interrupt_solve;
   PrimalDualHybridGradientParams params;
-  params.mutable_termination_criteria()->set_eps_optimal_absolute(0.0);
-  params.mutable_termination_criteria()->set_eps_optimal_relative(0.0);
+  params.mutable_termination_criteria()
+      ->mutable_simple_optimality_criteria()
+      ->set_eps_optimal_absolute(0.0);
+  params.mutable_termination_criteria()
+      ->mutable_simple_optimality_criteria()
+      ->set_eps_optimal_relative(0.0);
 
   interrupt_solve.store(false);
   auto callback = [&](const IterationCallbackInfo& info) {
@@ -1293,8 +1314,12 @@ TEST(PrimalDualHybridGradientTest, RespectsInterruptFromCallback) {
 TEST(PrimalDualHybridGradientTest, IgnoresFalseInterrupt) {
   std::atomic<bool> interrupt_solve;
   PrimalDualHybridGradientParams params;
-  params.mutable_termination_criteria()->set_eps_optimal_absolute(0.0);
-  params.mutable_termination_criteria()->set_eps_optimal_relative(0.0);
+  params.mutable_termination_criteria()
+      ->mutable_simple_optimality_criteria()
+      ->set_eps_optimal_absolute(0.0);
+  params.mutable_termination_criteria()
+      ->mutable_simple_optimality_criteria()
+      ->set_eps_optimal_relative(0.0);
   params.mutable_termination_criteria()->set_kkt_matrix_pass_limit(1);
 
   interrupt_solve.store(false);
@@ -1328,6 +1353,37 @@ TEST(PrimalDualHybridGradientTest, HugeNumShards) {
             TERMINATION_REASON_ITERATION_LIMIT);
 }
 
+TEST(PrimalDualHybridGradientTest, DetailedTerminationCriteria) {
+  const int iteration_upperbound = 300;
+  PrimalDualHybridGradientParams params = CreateSolverParams(
+      iteration_upperbound,
+      /*eps_optimal_absolute=*/1.0e-5, /*enable_scaling=*/true,
+      /*num_threads=*/4, /*use_iteration_limit=*/false,
+      /*use_malitsky_pock_linesearch=*/false,
+      /*use_diagonal_qp_trust_region_solver=*/false);
+  params.set_major_iteration_frequency(60);
+  params.mutable_termination_criteria()->clear_simple_optimality_criteria();
+  auto* opt_criteria = params.mutable_termination_criteria()
+                           ->mutable_detailed_optimality_criteria();
+  opt_criteria->set_eps_optimal_primal_residual_absolute(1e-5);
+  opt_criteria->set_eps_optimal_primal_residual_relative(0.0);
+  opt_criteria->set_eps_optimal_dual_residual_absolute(1e-5);
+  opt_criteria->set_eps_optimal_dual_residual_relative(0.0);
+  opt_criteria->set_eps_optimal_objective_gap_absolute(1e-5);
+  opt_criteria->set_eps_optimal_objective_gap_relative(0.0);
+
+  SolverResult output = PrimalDualHybridGradient(TinyLp(), params);
+  VerifyTerminationReasonAndIterationCount(params, output,
+                                           /*use_iteration_limit=*/false);
+  VerifyObjectiveValues(output, -1.0, 1.0e-4);
+  EXPECT_THAT(output.primal_solution,
+              EigenArrayNear<double>({1, 0, 6, 2}, 1.0e-4));
+  EXPECT_THAT(output.dual_solution,
+              EigenArrayNear<double>({0.5, 4.0, 0.0}, 1.0e-4));
+  EXPECT_THAT(output.reduced_costs,
+              EigenArrayNear<double>({0.0, 1.5, -3.5, 0.0}, 1.0e-4));
+}
+
 // Verifies that the primal and dual solution satisfy the bounds constraints.
 // This function uses ASSERTS rather than EXPECTs because it's used with large
 // QPs that would spam the logs otherwise.

ortools/pdlp/quadratic_program.cc

@@ -15,7 +15,6 @@
 
 #include <algorithm>
 #include <cstdint>
-#include <iterator>
 #include <limits>
 #include <string>
 #include <tuple>
@@ -27,7 +26,7 @@
 #include "absl/status/status.h"
 #include "absl/status/statusor.h"
 #include "absl/strings/str_cat.h"
-#include "ortools/base/logging.h"
+#include "ortools/base/check.h"
 #include "ortools/base/status_macros.h"
 #include "ortools/linear_solver/linear_solver.pb.h"
 

ortools/pdlp/quadratic_program.h

@@ -25,7 +25,6 @@
 #include "Eigen/SparseCore"
 #include "absl/status/status.h"
 #include "absl/status/statusor.h"
-#include "ortools/base/logging.h"
 #include "ortools/linear_solver/linear_solver.pb.h"
 
 namespace operations_research::pdlp {