ortools-clone/ortools/pdlp/primal_dual_hybrid_gradient_test.cc

// Copyright 2010-2025 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/pdlp/primal_dual_hybrid_gradient.h"

#include <atomic>
#include <cmath>
#include <cstdint>
#include <limits>
#include <optional>
#include <string>
#include <tuple>
#include <utility>
#include <vector>

#include "Eigen/Core"
#include "Eigen/SparseCore"
#include "absl/container/flat_hash_map.h"
#include "absl/log/check.h"
#include "absl/log/log.h"
#include "absl/status/statusor.h"
#include "absl/strings/str_cat.h"
#include "gtest/gtest.h"
#include "ortools/base/gmock.h"
#include "ortools/glop/parameters.pb.h"
#include "ortools/linear_solver/linear_solver.pb.h"
#include "ortools/lp_data/lp_data.h"
#include "ortools/lp_data/lp_types.h"
#include "ortools/pdlp/iteration_stats.h"
#include "ortools/pdlp/quadratic_program.h"
#include "ortools/pdlp/quadratic_program_io.h"
#include "ortools/pdlp/sharded_quadratic_program.h"
#include "ortools/pdlp/solve_log.pb.h"
#include "ortools/pdlp/solvers.pb.h"
#include "ortools/pdlp/termination.h"
#include "ortools/pdlp/test_util.h"

namespace operations_research::pdlp {
namespace {

using ::Eigen::VectorXd;
using ::operations_research::glop::ConstraintStatus;
using ::operations_research::glop::VariableStatus;
using ::testing::AnyOf;
using ::testing::DoubleNear;
using ::testing::ElementsAre;
using ::testing::Eq;
using ::testing::Ge;
using ::testing::HasSubstr;
using ::testing::IsEmpty;
using ::testing::Not;
using ::testing::Pair;
using ::testing::SizeIs;
using ::testing::UnorderedElementsAre;

const double kInfinity = std::numeric_limits<double>::infinity();

PrimalDualHybridGradientParams CreateSolverParams(
    const int iteration_limit, const double eps_optimal_absolute,
    const bool enable_scaling, const int num_threads,
    const bool use_iteration_limit, const bool use_malitsky_pock_linesearch,
    const bool use_diagonal_qp_trust_region_solver) {
  PrimalDualHybridGradientParams params;
  if (!enable_scaling) {
    params.set_l2_norm_rescaling(false);
    params.set_l_inf_ruiz_iterations(0);
  }
  if (use_malitsky_pock_linesearch) {
    params.set_linesearch_rule(
        PrimalDualHybridGradientParams::MALITSKY_POCK_LINESEARCH_RULE);
  }

  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()
        ->mutable_simple_optimality_criteria()
        ->set_eps_optimal_absolute(0.0);
  } else {
    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);
    params.set_diagonal_qp_trust_region_solver_tolerance(1.0e-8);
  }
  params.set_num_threads(num_threads);
  // Protect against infinite loops if the termination criteria fail.
  params.mutable_termination_criteria()->set_kkt_matrix_pass_limit(1'000'000.0);
  return params;
}

// Verifies expected termination reason and iteration count for an instance
// where an optimal solution exists.
// `params` must have been generated by `CreateSolverParams()` with the same
// `use_iteration_limit`.
// Intended usage:
//   const bool use_iteration_limit = ...;
//   PrimalDualHybridGradientParams params =
//       CreateSolverParams(..., use_iteration_limit, ...);
//    SolverResult output = PrimalDualHybridGradient(..., params);
//    VerifyTerminationReasonAndIterationCount(params, output,
//                                             use_iteration_limit);
void VerifyTerminationReasonAndIterationCount(
    const PrimalDualHybridGradientParams& params, const SolverResult& output,
    const bool use_iteration_limit) {
  if (use_iteration_limit) {
    // When a PDHG step has zero length PDHG can no longer make progress and
    // hence terminates immediately. In theory a zero step length implies
    // optimality but in practice PDHG terminates with a reason of OPTIMAL if
    // the optimality checks pass and NUMERICAL_ERROR otherwise.
    // When `use_iteration_limit==true`, `CreateSolverParams()` sets all the
    // epsilons to 0, which makes the optimality checks harder to pass but not
    // impossible. Both OPTIMAL and NUMERICAL_ERROR are therefore ok termination
    // reasons.
    EXPECT_THAT(
        output.solve_log.termination_reason(),
        AnyOf(TERMINATION_REASON_ITERATION_LIMIT,
              TERMINATION_REASON_NUMERICAL_ERROR, TERMINATION_REASON_OPTIMAL));
    if (output.solve_log.termination_reason() ==
        TERMINATION_REASON_ITERATION_LIMIT) {
      EXPECT_EQ(output.solve_log.iteration_count(),
                params.termination_criteria().iteration_limit());
    } else {
      EXPECT_LE(output.solve_log.iteration_count(),
                params.termination_criteria().iteration_limit());
    }
  } else {
    EXPECT_EQ(output.solve_log.termination_reason(),
              TERMINATION_REASON_OPTIMAL);
    EXPECT_LE(output.solve_log.iteration_count(),
              params.termination_criteria().iteration_limit());
  }
}

// Verifies the primal and dual objective values.
void VerifyObjectiveValues(const SolverResult& result,
                           const double objective_value,
                           const double tolerance) {
  const auto& convergence_info = GetConvergenceInformation(
      result.solve_log.solution_stats(), result.solve_log.solution_type());
  ASSERT_TRUE(convergence_info.has_value());
  EXPECT_THAT(convergence_info->primal_objective(),
              DoubleNear(objective_value, tolerance));
  EXPECT_THAT(convergence_info->dual_objective(),
              DoubleNear(objective_value, tolerance));
}

class PrimalDualHybridGradientLPTest
    : public testing::TestWithParam<
          std::tuple</*enable_scaling=*/bool, /*num_threads=*/int,
                     /*use_iteration_limit=*/bool,
                     /*use_malitsky_pock_linesearch=*/bool>> {
 protected:
  PrimalDualHybridGradientParams CreateSolverParamsForFixture(
      const int iteration_limit, const double eps_optimal_absolute) {
    const auto [enable_scaling, num_threads, use_iteration_limit,
                use_malitsky_pock_linesearch] = GetParam();
    return CreateSolverParams(iteration_limit, eps_optimal_absolute,
                              enable_scaling, num_threads, use_iteration_limit,
                              use_malitsky_pock_linesearch,
                              /*use_diagonal_qp_trust_region_solver=*/false);
  }

  void VerifyTerminationReasonAndIterationCountForFixture(
      const PrimalDualHybridGradientParams& params,
      const SolverResult& output) {
    const auto [enable_scaling, num_threads, use_iteration_limit,
                use_malitsky_pock_linesearch] = GetParam();
    VerifyTerminationReasonAndIterationCount(params, output,
                                             use_iteration_limit);
  }
};

class PrimalDualHybridGradientDiagonalQPTest
    : public testing::TestWithParam<
          std::tuple</*enable_scaling=*/bool, /*num_threads=*/int,
                     /*use_iteration_limit=*/bool,
                     /*use_malitsky_pock_linesearch=*/bool,
                     /*use_diagonal_qp_trust_region_solver=*/bool>> {
 protected:
  PrimalDualHybridGradientParams CreateSolverParamsForFixture(
      const int iteration_limit, const double eps_optimal_absolute) {
    const auto [enable_scaling, num_threads, use_iteration_limit,
                use_malitsky_pock_linesearch,
                use_diagonal_qp_trust_region_solver] = GetParam();
    return CreateSolverParams(iteration_limit, eps_optimal_absolute,
                              enable_scaling, num_threads, use_iteration_limit,
                              use_malitsky_pock_linesearch,
                              use_diagonal_qp_trust_region_solver);
  }
  void VerifyTerminationReasonAndIterationCountForFixture(
      const PrimalDualHybridGradientParams& params,
      const SolverResult& output) {
    const auto [enable_scaling, num_threads, use_iteration_limit,
                use_malitsky_pock_linesearch,
                use_diagonal_qp_trust_region_solver] = GetParam();
    VerifyTerminationReasonAndIterationCount(params, output,
                                             use_iteration_limit);
  }
};

class PrimalDualHybridGradientVerbosityTest
    : public testing::TestWithParam</*verbosity_level=*/int> {};

class PresolveDualScalingTest
    : public testing::TestWithParam<
          std::tuple</*Dualize=*/bool,
                     /*NegateAndScaleObjective=*/bool>> {};

INSTANTIATE_TEST_SUITE_P(
    QP, PrimalDualHybridGradientDiagonalQPTest,
    testing::Combine(testing::Bool(), testing::Values(1, 4), testing::Bool(),
                     testing::Bool(), testing::Bool()),
    [](const testing::TestParamInfo<
        PrimalDualHybridGradientDiagonalQPTest::ParamType>& info) {
      return absl::StrCat(
          std::get<0>(info.param) ? "Scaling" : "NoScaling", "_",
          std::get<1>(info.param), "Threads_",
          std::get<2>(info.param) ? "IterationLimit" : "NoIterationLimit", "_",
          std::get<3>(info.param) ? "MalitskyPockLinesearch"
                                  : "AdaptiveLinesearch",
          "_", std::get<4>(info.param) ? "TRSolverDiag" : "TRSolverNoDiag");
    });

INSTANTIATE_TEST_SUITE_P(
    LP, PrimalDualHybridGradientLPTest,
    testing::Combine(testing::Bool(), testing::Values(1, 4), testing::Bool(),
                     testing::Bool()),
    [](const testing::TestParamInfo<PrimalDualHybridGradientLPTest::ParamType>&
           info) {
      return absl::StrCat(
          std::get<0>(info.param) ? "Scaling" : "NoScaling", "_",
          std::get<1>(info.param), "Threads_",
          std::get<2>(info.param) ? "IterationLimit" : "NoIterationLimit", "_",
          std::get<3>(info.param) ? "MalitskyPockLinesearch"
                                  : "AdaptiveLinesearch");
    });

INSTANTIATE_TEST_SUITE_P(Verbosity, PrimalDualHybridGradientVerbosityTest,
                         testing::Values(0, 1, 2, 3, 4));

INSTANTIATE_TEST_SUITE_P(
    PresolveDualScaling, PresolveDualScalingTest,
    testing::Combine(testing::Bool(), testing::Bool()),
    [](const testing::TestParamInfo<PresolveDualScalingTest::ParamType>& info) {
      return absl::StrCat(std::get<1>(info.param) ? "Dualize" : "NoDualize",
                          std::get<0>(info.param) ? "NegateAndScaleObjective"
                                                  : "NoObjectiveScaling");
    });

TEST_P(PrimalDualHybridGradientLPTest, UnboundedVariables) {
  const int iteration_upperbound = 980;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-7);
  params.set_major_iteration_frequency(100);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  VerifyTerminationReasonAndIterationCountForFixture(params, output);
  VerifyObjectiveValues(output, -34.0, 1.0e-6);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({-1, 8, 1, 2.5}, 1.0e-4));
  EXPECT_THAT(output.dual_solution,
              EigenArrayNear<double>({-2, 0, 2.375, 2.0 / 3}, 1.0e-4));

  EXPECT_EQ(output.solve_log.original_problem_stats().num_variables(), 4);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_variables(), 4);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_constraints(), 4);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_constraints(), 4);
}

TEST_P(PrimalDualHybridGradientLPTest, Tiny) {
  const int iteration_upperbound = 300;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-5);
  params.set_major_iteration_frequency(60);

  SolverResult output = PrimalDualHybridGradient(TinyLp(), params);
  VerifyTerminationReasonAndIterationCountForFixture(params, output);
  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));
  EXPECT_EQ(output.solve_log.original_problem_stats().num_variables(), 4);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_variables(), 4);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_constraints(), 3);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_constraints(), 3);
}

TEST_P(PrimalDualHybridGradientLPTest, CorrelationClusteringOne) {
  const int iteration_upperbound = 9;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-10);
  params.set_major_iteration_frequency(2);

  SolverResult output =
      PrimalDualHybridGradient(CorrelationClusteringLp(), params);
  VerifyTerminationReasonAndIterationCountForFixture(params, output);
  VerifyObjectiveValues(output, 1.0, 1.0e-14);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({1, 1, 0, 1, 0, 0}, 1.0e-14));
  ASSERT_EQ(output.dual_solution.size(), 3);
  // There are multiple optimal dual solutions.
  EXPECT_GE(output.dual_solution[0], 0);
  EXPECT_GE(output.dual_solution[1], 0);
  EXPECT_GE(output.dual_solution[2], 0);
  EXPECT_GE(output.dual_solution[0] + output.dual_solution[1], 1 - 1.0e-14);

  const auto& convergence_information = GetConvergenceInformation(
      output.solve_log.solution_stats(), output.solve_log.solution_type());
  ASSERT_TRUE(convergence_information.has_value());
  EXPECT_THAT(convergence_information->corrected_dual_objective(),
              DoubleNear(1, 1.0e-14));
  EXPECT_EQ(output.solve_log.original_problem_stats().num_variables(), 6);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_variables(), 6);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_constraints(), 3);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_constraints(), 3);
}

TEST_P(PrimalDualHybridGradientLPTest, CorrelationClusteringStar) {
  const int iteration_upperbound = 45;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);
  params.set_major_iteration_frequency(5);

  SolverResult output =
      PrimalDualHybridGradient(CorrelationClusteringStarLp(), params);
  VerifyTerminationReasonAndIterationCountForFixture(params, output);
  VerifyObjectiveValues(output, 1.5, 1.0e-6);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({0.5, 0.5, 0.5, 0.0, 0.0, 0.0}, 1.0e-6));
  EXPECT_THAT(output.dual_solution,
              EigenArrayNear<double>({0.5, 0.5, 0.5}, 1.0e-6));
  EXPECT_EQ(output.solve_log.original_problem_stats().num_variables(), 6);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_variables(), 6);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_constraints(), 3);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_constraints(), 3);
}

// A double-sided constraint l <= a^T x <= u where neither constraint is tight
// at optimum could cause the trust region solver to malfunction if we picked
// the wrong dual subgradient. This test verifies that we solve an instance with
// such a constraint quickly to high accuracy.
TEST_P(PrimalDualHybridGradientLPTest, InactiveTwoSidedConstraint) {
  const int iteration_upperbound = 500;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-8);
  params.set_major_iteration_frequency(60);

  QuadraticProgram qp = TestLp();
  // This makes this constraint double-sided and inactive at the optimal
  // solution.
  qp.constraint_lower_bounds[1] = -10;
  SolverResult output = PrimalDualHybridGradient(qp, params);
  VerifyTerminationReasonAndIterationCountForFixture(params, output);
  VerifyObjectiveValues(output, -34.0, 1.0e-6);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({-1, 8, 1, 2.5}, 1.0e-7));
  EXPECT_THAT(output.dual_solution,
              EigenArrayNear<double>({-2.0, 0.0, 2.375, 2.0 / 3}, 1.0e-7));
}

TEST_P(PrimalDualHybridGradientLPTest, InfeasiblePrimal) {
  // This value for `iteration_upperbound` is particularly necessary for
  // Malistsky and Pock. The adaptive rule detects infeasibility in less than
  // 500 iterations.
  const int iteration_upperbound = 2000;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);
  params.set_major_iteration_frequency(5);
  params.mutable_termination_criteria()->set_eps_primal_infeasible(1.0e-6);

  SolverResult output =
      PrimalDualHybridGradient(SmallPrimalInfeasibleLp(), params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_PRIMAL_INFEASIBLE);
  const auto& dual = output.dual_solution;
  // The following two conditions check if the certificate is correct. For this
  // problem the set of infeasibility certificates is equal to all the rays of
  // the form -alpha * (1, 1) with alpha positive.
  EXPECT_THAT(dual[0] / dual[1], DoubleNear(1, 1.0e-6));
  EXPECT_LT(dual[1], 0.0);
  // The reduced costs should be approximately zero. However, a small relative
  // difference between `dual[0]` and `dual[1]` could translate to a large
  // absolute difference, and hence large reduced costs. The following test uses
  // the exact formula to make sure we're not adding the objective vector to the
  // reduced costs.
  EXPECT_THAT(
      output.reduced_costs,
      EigenArrayNear<double>({dual[1] - dual[0], dual[0] - dual[1]}, 1.0e-6));
  EXPECT_LE(output.solve_log.iteration_count(), iteration_upperbound);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_variables(), 2);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_variables(), 2);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_constraints(), 2);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_constraints(), 2);
}

TEST_P(PrimalDualHybridGradientLPTest, InfeasibleDual) {
  const int iteration_upperbound = 500;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);
  params.set_major_iteration_frequency(5);

  SolverResult output =
      PrimalDualHybridGradient(SmallDualInfeasibleLp(), params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_DUAL_INFEASIBLE);
  // The following two conditions check if the certificate is correct. For this
  // problem the set of infeasibility certificates is equal to all the rays of
  // the form alpha * (1, 1) with alpha positive.
  EXPECT_THAT(output.primal_solution[0] / output.primal_solution[1],
              DoubleNear(1, 1.0e-6));
  EXPECT_GT(output.primal_solution[1], 0.0);
  EXPECT_LE(output.solve_log.iteration_count(), iteration_upperbound);
}

TEST_P(PrimalDualHybridGradientLPTest, InfeasiblePrimalWithReducedCosts) {
  // A trivial LP that has reduced costs within the dual ray.
  //     min x
  //     Constraint: 2 <= x
  //     Variable: 0 <= x <= 1
  QuadraticProgram lp(1, 1);
  lp.objective_vector = Eigen::VectorXd{{1}};
  lp.constraint_lower_bounds = Eigen::VectorXd{{2}};
  lp.constraint_upper_bounds =
      Eigen::VectorXd{{std::numeric_limits<double>::infinity()}};
  lp.variable_lower_bounds = Eigen::VectorXd{{0}};
  lp.variable_upper_bounds = Eigen::VectorXd{{1}};
  lp.constraint_matrix.coeffRef(0, 0) = 1.0;
  lp.constraint_matrix.makeCompressed();
  const int iteration_upperbound = 100;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);
  params.set_major_iteration_frequency(5);
  params.mutable_termination_criteria()->set_eps_primal_infeasible(1.0e-6);

  SolverResult output = PrimalDualHybridGradient(lp, params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_PRIMAL_INFEASIBLE);
  const auto& dual = output.dual_solution;
  EXPECT_GT(dual[0], 0.0);
  EXPECT_EQ(output.reduced_costs[0], -dual[0]);
  EXPECT_LT(output.solve_log.iteration_count(), iteration_upperbound);
}

TEST_P(PrimalDualHybridGradientLPTest, InfeasiblePrimalDual) {
  const int iteration_upperbound = 600;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);
  // Adaptive restarts are disabled because they unexpectedly perform worse on
  // this instance.
  params.set_restart_strategy(PrimalDualHybridGradientParams::NO_RESTARTS);
  params.set_major_iteration_frequency(5);

  SolverResult output =
      PrimalDualHybridGradient(SmallPrimalDualInfeasibleLp(), params);
  EXPECT_THAT(output.solve_log.termination_reason(),
              AnyOf(Eq(TERMINATION_REASON_DUAL_INFEASIBLE),
                    Eq(TERMINATION_REASON_PRIMAL_INFEASIBLE)));
  EXPECT_LE(output.solve_log.iteration_count(), iteration_upperbound);
}

TEST_P(PrimalDualHybridGradientDiagonalQPTest, DiagonalQp1) {
  const int iteration_upperbound = 96;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);
  params.set_major_iteration_frequency(12);

  SolverResult output = PrimalDualHybridGradient(TestDiagonalQp1(), params);
  VerifyTerminationReasonAndIterationCountForFixture(params, output);
  VerifyObjectiveValues(output, 6.0, 1.0e-6);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({1.0, 0.0}, 1.0e-6));
  EXPECT_THAT(output.dual_solution,
              EigenArrayNear(Eigen::ArrayXd::Constant(1, -1.0), 1.0e-6));
  EXPECT_THAT(output.reduced_costs, EigenArrayNear<double>({4.0, 0.0}, 1.0e-6));
  EXPECT_EQ(output.solve_log.original_problem_stats().num_variables(), 2);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_variables(), 2);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_constraints(), 1);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_constraints(), 1);
}

TEST_P(PrimalDualHybridGradientDiagonalQPTest, DiagonalQp2) {
  const int iteration_upperbound = 240;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);
  params.set_major_iteration_frequency(12);

  SolverResult output = PrimalDualHybridGradient(TestDiagonalQp2(), params);
  VerifyTerminationReasonAndIterationCountForFixture(params, output);
  VerifyObjectiveValues(output, -5.0, 1.0e-6);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({3.0, 1.0}, 1.0e-6));
  EXPECT_THAT(output.dual_solution, ElementsAre(DoubleNear(0.0, 1.0e-6)));
  EXPECT_THAT(output.reduced_costs, EigenArrayNear<double>({0.0, 0.0}, 1.0e-6));
}

TEST_P(PrimalDualHybridGradientDiagonalQPTest, DiagonalQp3) {
  const int iteration_upperbound = 300;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);
  params.set_major_iteration_frequency(15);

  SolverResult output = PrimalDualHybridGradient(TestDiagonalQp3(), params);
  VerifyTerminationReasonAndIterationCountForFixture(params, output);
  VerifyObjectiveValues(output, 2.0, 1.0e-6);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({2.0, 0.0, 1.0}, 1.0e-6));
  EXPECT_THAT(output.dual_solution,
              EigenArrayNear<double>({-1.0, 1.0}, 1.0e-6));
  EXPECT_THAT(output.reduced_costs, EigenArrayNear<double>({0, 0, 0}, 1.0e-6));
}

// This is like `DiagonalQp1` except it starts with a near-optimal solution and
// uses a shorter iteration limit.
TEST_P(PrimalDualHybridGradientDiagonalQPTest, QpWarmStart) {
  const int iteration_upperbound = 35;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);
  params.set_major_iteration_frequency(5);
  // Disable primal weight updating. In a warm-start situation, the primal
  // weight should be carried over. In this test, the initial primal weight of 1
  // is reasonable because the distance from the starting point to the primal
  // and dual optimal solutions are about the same.
  params.set_primal_weight_update_smoothing(0.0);

  const PrimalAndDualSolution initial_solution = {
      .primal_solution = Eigen::VectorXd{{0.999, 0.001}},
      .dual_solution = Eigen::VectorXd{{-0.999}},
  };
  SolverResult output = PrimalDualHybridGradient(TestDiagonalQp1(), params,
                                                 std::move(initial_solution));
  VerifyTerminationReasonAndIterationCountForFixture(params, output);
  VerifyObjectiveValues(output, 6.0, 1.0e-6);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({1.0, 0.0}, 1.0e-6));
  EXPECT_THAT(output.dual_solution,
              EigenArrayNear(Eigen::ArrayXd::Constant(1, -1.0), 1.0e-6));
  EXPECT_THAT(output.reduced_costs, EigenArrayNear<double>({4.0, 0.0}, 1.0e-6));
}

// Tests an LP with no constraints.
TEST_P(PrimalDualHybridGradientLPTest, LpWithoutConstraints) {
  const int iteration_upperbound = 2;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);

  QuadraticProgram qp(3, 0);
  qp.variable_lower_bounds = Eigen::VectorXd{{-1, -kInfinity, -2}};
  qp.variable_upper_bounds = Eigen::VectorXd{{kInfinity, 4, 10}};
  qp.objective_vector = Eigen::VectorXd{{1, -1, 2}};

  SolverResult output = PrimalDualHybridGradient(qp, params);
  VerifyTerminationReasonAndIterationCountForFixture(params, output);
  VerifyObjectiveValues(output, -9.0, 1.0e-6);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({-1, 4, -2}, 1.0e-6));
  EXPECT_THAT(output.dual_solution, SizeIs(0));
  EXPECT_EQ(output.solve_log.original_problem_stats().num_variables(), 3);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_variables(), 3);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_constraints(), 0);
  EXPECT_EQ(output.solve_log.preprocessed_problem_stats().num_constraints(), 0);
}

// Tests an LP with no variables.
TEST_P(PrimalDualHybridGradientLPTest, LpWithoutVariables) {
  const int iteration_upperbound = 2;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);

  QuadraticProgram qp(0, 3);
  qp.constraint_lower_bounds = Eigen::VectorXd{{-1, -kInfinity, -2}};
  qp.constraint_upper_bounds = Eigen::VectorXd{{kInfinity, 4, 10}};

  SolverResult output = PrimalDualHybridGradient(qp, params);
  VerifyTerminationReasonAndIterationCountForFixture(params, output);
  VerifyObjectiveValues(output, 0.0, 1.0e-6);
  EXPECT_THAT(output.primal_solution, SizeIs(0));
  EXPECT_THAT(output.dual_solution, EigenArrayNear<double>({0, 0, 0}, 1.0e-6));
  EXPECT_EQ(output.solve_log.original_problem_stats().num_variables(), 0);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_variables(), 0);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_constraints(), 3);
  EXPECT_EQ(output.solve_log.preprocessed_problem_stats().num_constraints(), 3);
}

TEST_P(PrimalDualHybridGradientLPTest, LpWithOnlyFixedVariable) {
  const int iteration_upperbound = 2;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/0.0);

  QuadraticProgram qp(1, 0);
  qp.variable_lower_bounds = Eigen::VectorXd{{1}};
  qp.variable_upper_bounds = Eigen::VectorXd{{1}};
  qp.objective_vector = Eigen::VectorXd{{1}};

  SolverResult output = PrimalDualHybridGradient(qp, params);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear(Eigen::ArrayXd::Constant(1, 1.0), 1.0e-6));
  EXPECT_THAT(output.dual_solution, testing::SizeIs(0));
  EXPECT_LE(output.solve_log.iteration_count(), iteration_upperbound);
}

TEST_P(PrimalDualHybridGradientLPTest, InfeasibleLpWithoutVariables) {
  const int iteration_upperbound = 2;
  PrimalDualHybridGradientParams params =
      CreateSolverParamsForFixture(iteration_upperbound,
                                   /*eps_optimal_absolute=*/1.0e-6);

  QuadraticProgram qp(0, 1);
  qp.constraint_lower_bounds = Eigen::VectorXd{{-1}};
  qp.constraint_upper_bounds = Eigen::VectorXd{{-1}};

  SolverResult output = PrimalDualHybridGradient(qp, params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_PRIMAL_INFEASIBLE);
  EXPECT_THAT(output.primal_solution, SizeIs(0));
  EXPECT_LT(output.dual_solution[0], 0.0);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_variables(), 0);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_variables(), 0);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_constraints(), 1);
  EXPECT_EQ(output.solve_log.preprocessed_problem_stats().num_constraints(), 1);
}

PrimalDualHybridGradientParams ParamsWithNoLimits() {
  PrimalDualHybridGradientParams params;
  // This disables the termination limits. A termination criteria must be set
  // for the solver to terminate.
  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;
}

TEST(PrimalDualHybridGradientTest, ClearsRunningAverage) {
  // An arbitrarily chosen major iteration frequency.
  const int major_iteration_frequency = 17;
  const int iteration_limit = 100;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_termination_check_frequency(1);
  params.set_major_iteration_frequency(major_iteration_frequency);
  params.set_restart_strategy(PrimalDualHybridGradientParams::NO_RESTARTS);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  ASSERT_EQ(output.solve_log.iteration_count(), iteration_limit);
  // The first entry in `iteration_stats` corresponds to the starting point.
  ASSERT_EQ(output.solve_log.iteration_stats_size(), iteration_limit + 1);

  for (int i = 0; i < output.solve_log.iteration_stats_size(); ++i) {
    const auto& stats = output.solve_log.iteration_stats(i);
    const int iterations_completed = stats.iteration_number();
    EXPECT_EQ(iterations_completed, i);
    if (iterations_completed == 0) {
      EXPECT_EQ(stats.restart_used(), RESTART_CHOICE_NO_RESTART);
    } else if (iterations_completed % major_iteration_frequency == 0) {
      EXPECT_EQ(stats.restart_used(), RESTART_CHOICE_WEIGHTED_AVERAGE_RESET)
          << "iteration = " << i;
    } else {
      EXPECT_EQ(stats.restart_used(), RESTART_CHOICE_NO_RESTART)
          << "iteration = " << i;
    }
  }
}

TEST(PrimalDualHybridGradientTest, RestartsEveryMajorIteration) {
  // An arbitrarily chosen major iteration frequency.
  const int major_iteration_frequency = 17;
  const int iteration_limit = 100;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_termination_check_frequency(1);
  params.set_major_iteration_frequency(major_iteration_frequency);
  params.set_restart_strategy(
      PrimalDualHybridGradientParams::EVERY_MAJOR_ITERATION);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  ASSERT_EQ(output.solve_log.iteration_count(), iteration_limit);
  // The first entry in `iteration_stats` corresponds to the starting point.
  ASSERT_EQ(output.solve_log.iteration_stats_size(), iteration_limit + 1);

  for (int i = 0; i < output.solve_log.iteration_stats_size(); ++i) {
    const auto& stats = output.solve_log.iteration_stats(i);
    const int iterations_completed = stats.iteration_number();
    EXPECT_EQ(iterations_completed, i);
    if (iterations_completed == 0) {
      EXPECT_EQ(stats.restart_used(), RESTART_CHOICE_NO_RESTART);
    } else if (iterations_completed % major_iteration_frequency == 0) {
      EXPECT_EQ(stats.restart_used(), RESTART_CHOICE_RESTART_TO_AVERAGE)
          << "iteration = " << i;
    } else {
      EXPECT_EQ(stats.restart_used(), RESTART_CHOICE_NO_RESTART)
          << "iteration = " << i;
    }
  }
}

TEST(PrimalDualHybridGradientTest, SolveLogIncludesNameForNamedQP) {
  PrimalDualHybridGradientParams params;
  params.mutable_termination_criteria()->set_iteration_limit(1);

  QuadraticProgram test_lp = TestLp();
  test_lp.problem_name = "Test LP";

  SolverResult output = PrimalDualHybridGradient(test_lp, params);
  EXPECT_EQ(output.solve_log.instance_name(), "Test LP");
}

TEST(PrimalDualHybridGradientTest, SolveLogOmitsNameForUnnamedQP) {
  PrimalDualHybridGradientParams params;
  params.mutable_termination_criteria()->set_iteration_limit(1);

  QuadraticProgram unnamed_test_lp = TestLp();

  SolverResult output = PrimalDualHybridGradient(unnamed_test_lp, params);
  EXPECT_FALSE(output.solve_log.has_instance_name());
}

TEST(PrimalDualHybridGradientTest, SolveLogIncludesParameters) {
  PrimalDualHybridGradientParams params;
  params.mutable_termination_criteria()->set_iteration_limit(1);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  EXPECT_EQ(output.solve_log.params().termination_criteria().iteration_limit(),
            1);
}

TEST(PrimalDualHybridGradientTest, AdaptiveDistanceBasedRestartsWorkOnTestLp) {
  PrimalDualHybridGradientParams params;
  params.set_major_iteration_frequency(16);
  params.mutable_termination_criteria()->set_iteration_limit(128);
  params.set_restart_strategy(
      PrimalDualHybridGradientParams::ADAPTIVE_DISTANCE_BASED);
  // Low restart threshold.
  params.set_necessary_reduction_for_restart(0.99);
  SolverResult output = PrimalDualHybridGradient(TestLp(), params);

  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({-1, 8, 1, 2.5}, 1.0e-4));
  EXPECT_THAT(output.dual_solution,
              EigenArrayNear<double>({-2, 0, 2.375, 2.0 / 3}, 1.0e-4));
  const auto convergence_info = GetConvergenceInformation(
      output.solve_log.solution_stats(), output.solve_log.solution_type());
  ASSERT_TRUE(convergence_info.has_value());
  EXPECT_THAT(convergence_info->primal_objective(), DoubleNear(-34.0, 1.0e-4));
  EXPECT_THAT(convergence_info->dual_objective(), DoubleNear(-34.0, 1.0e-4));
}

TEST(PrimalDualHybridGradientTest, AdaptiveDistanceBasedRestartsWorkOnTestQp) {
  PrimalDualHybridGradientParams params;
  params.set_major_iteration_frequency(16);
  params.mutable_termination_criteria()->set_iteration_limit(128);
  params.set_restart_strategy(
      PrimalDualHybridGradientParams::ADAPTIVE_DISTANCE_BASED);
  // Low restart threshold.
  params.set_necessary_reduction_for_restart(0.99);
  SolverResult output = PrimalDualHybridGradient(TestDiagonalQp1(), params);

  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST(PrimalDualHybridGradientTest, AdaptiveDistanceBasedRestartsToAverage) {
  // An arbitrarily chosen major iteration frequency.
  const int major_iteration_frequency = 13;
  const int iteration_limit = 100;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_termination_check_frequency(1);
  params.set_major_iteration_frequency(major_iteration_frequency);
  params.set_restart_strategy(
      PrimalDualHybridGradientParams::ADAPTIVE_DISTANCE_BASED);
  params.set_necessary_reduction_for_restart(0.75);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  ASSERT_EQ(output.solve_log.iteration_count(), iteration_limit);
  // The first entry in `iteration_stats` corresponds to the starting point.
  ASSERT_EQ(output.solve_log.iteration_stats_size(), iteration_limit + 1);

  for (int i = 0; i < output.solve_log.iteration_stats_size(); ++i) {
    const auto& stats = output.solve_log.iteration_stats(i);
    const int iterations_completed = stats.iteration_number();
    EXPECT_EQ(iterations_completed, i);
    if (iterations_completed == 0) {
      EXPECT_EQ(stats.restart_used(), RESTART_CHOICE_NO_RESTART);
    } else if (iterations_completed == major_iteration_frequency) {
      // An explicit restart should be triggered at the end of the first major
      // iteration.
      EXPECT_THAT(stats.restart_used(),
                  AnyOf(RESTART_CHOICE_RESTART_TO_AVERAGE,
                        RESTART_CHOICE_WEIGHTED_AVERAGE_RESET))
          << "iteration = " << i;
    } else if (iterations_completed % major_iteration_frequency != 0) {
      // No restarts should happen outside major iterations.
      EXPECT_EQ(stats.restart_used(), RESTART_CHOICE_NO_RESTART);
    }
  }
}

TEST(PrimalDualHybridGradientTest, PrimalWeightFrozen) {
  // An arbitrarily chosen major iteration frequency.
  const int major_iteration_frequency = 17;
  const int iteration_limit = 100;
  const double initial_primal_weight = 1.5;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_major_iteration_frequency(major_iteration_frequency);
  params.set_restart_strategy(
      PrimalDualHybridGradientParams::EVERY_MAJOR_ITERATION);
  params.set_initial_primal_weight(initial_primal_weight);
  params.set_primal_weight_update_smoothing(0.0);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);

  for (const auto& stats : output.solve_log.iteration_stats()) {
    EXPECT_EQ(stats.primal_weight(), initial_primal_weight)
        << "iteration = " << stats.iteration_number();
  }
}

TEST(PrimalDualHybridGradientTest, ConstantStepSize) {
  const int iteration_limit = 100;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_termination_check_frequency(1);
  params.set_linesearch_rule(
      PrimalDualHybridGradientParams::CONSTANT_STEP_SIZE_RULE);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);

  ASSERT_FALSE(output.solve_log.iteration_stats().empty());
  const double initial_step_size =
      output.solve_log.iteration_stats(0).step_size();
  for (const auto& stats : output.solve_log.iteration_stats()) {
    EXPECT_EQ(stats.step_size(), initial_step_size)
        << "iteration = " << stats.iteration_number();
  }
}

TEST(PrimalDualHybridGradientTest, StepSizeScaling) {
  const int iteration_limit = 1;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_termination_check_frequency(1);
  params.set_linesearch_rule(
      PrimalDualHybridGradientParams::CONSTANT_STEP_SIZE_RULE);

  SolverResult unscaled_output = PrimalDualHybridGradient(TestLp(), params);

  ASSERT_FALSE(unscaled_output.solve_log.iteration_stats().empty());
  const double initial_step_size =
      unscaled_output.solve_log.iteration_stats(0).step_size();

  const double kStepSizeScaling = 0.5;
  params.set_initial_step_size_scaling(kStepSizeScaling);
  SolverResult scaled_output = PrimalDualHybridGradient(TestLp(), params);

  ASSERT_FALSE(scaled_output.solve_log.iteration_stats().empty());
  EXPECT_EQ(scaled_output.solve_log.iteration_stats(0).step_size(),
            initial_step_size * kStepSizeScaling);
}

// This verifies that `kkt_matrix_pass_limit` is checked every iteration.
TEST(PrimalDualHybridGradientTest, KktMatrixPassTermination) {
  const int kkt_matrix_pass_limit = 13;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_kkt_matrix_pass_limit(
      kkt_matrix_pass_limit);
  params.set_linesearch_rule(
      PrimalDualHybridGradientParams::CONSTANT_STEP_SIZE_RULE);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);

  ASSERT_FALSE(output.solve_log.iteration_stats().empty());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT);
  EXPECT_EQ(output.solve_log.solution_stats().cumulative_kkt_matrix_passes(),
            kkt_matrix_pass_limit);
}

TEST(PrimalDualHybridGradientTest,
     StatsAtEachIterationWithRecordIterationStatsOn) {
  // An arbitrarily chosen major iteration frequency.
  const int major_iteration_frequency = 17;
  const int iteration_limit = 100;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.set_record_iteration_stats(true);
  // This is required for `ConvergenceInformation`, `InfeasibilityInformation`,
  // and `PointMetadata` to be generated on each iteration.
  params.set_termination_check_frequency(1);
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_major_iteration_frequency(major_iteration_frequency);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  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),
              std::nullopt);
    EXPECT_NE(GetInfeasibilityInformation(stats, POINT_TYPE_CURRENT_ITERATE),
              std::nullopt);
    EXPECT_NE(GetPointMetadata(stats, POINT_TYPE_CURRENT_ITERATE),
              std::nullopt);
    if (stats.iteration_number() > 0) {
      EXPECT_NE(GetConvergenceInformation(stats, POINT_TYPE_AVERAGE_ITERATE),
                std::nullopt);
      EXPECT_NE(GetInfeasibilityInformation(stats, POINT_TYPE_AVERAGE_ITERATE),
                std::nullopt);
      EXPECT_NE(GetPointMetadata(stats, POINT_TYPE_AVERAGE_ITERATE),
                std::nullopt);

      EXPECT_NE(
          GetInfeasibilityInformation(stats, POINT_TYPE_ITERATE_DIFFERENCE),
          std::nullopt);
      EXPECT_NE(GetPointMetadata(stats, POINT_TYPE_ITERATE_DIFFERENCE),
                std::nullopt);
    }
  }
}

TEST(PrimalDualHybridGradientTest,
     NoIterationStatsWithRecordIterationStatsOff) {
  // An arbitrarily chosen major iteration frequency.
  const int major_iteration_frequency = 17;
  const int iteration_limit = 100;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.set_record_iteration_stats(false);
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_major_iteration_frequency(major_iteration_frequency);
  // Random projection seeds should have no effect when `record_iteration_stats`
  // is false.
  params.add_random_projection_seeds(1);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  EXPECT_EQ(output.solve_log.iteration_stats().size(), 0);
}

TEST(PrimalDualHybridGradientTest, NoRandomProjectionsIfNotRequested) {
  // An arbitrarily chosen major iteration frequency.
  const int major_iteration_frequency = 17;
  const int iteration_limit = 100;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.set_record_iteration_stats(true);
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_major_iteration_frequency(major_iteration_frequency);
  params.set_termination_check_frequency(1);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  EXPECT_EQ(output.solve_log.iteration_stats().size(), iteration_limit + 1);
  for (const auto& stats : output.solve_log.iteration_stats()) {
    EXPECT_THAT(stats.point_metadata(), Not(IsEmpty()));
    for (const auto& metadata : stats.point_metadata()) {
      EXPECT_THAT(metadata.random_primal_projections(), IsEmpty());
      EXPECT_THAT(metadata.random_dual_projections(), IsEmpty());
    }
  }
}

TEST(PrimalDualHybridGradientTest, HasRandomProjectionsIfRequested) {
  // An arbitrarily chosen major iteration frequency.
  const int major_iteration_frequency = 17;
  const int iteration_limit = 100;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.set_record_iteration_stats(true);
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_major_iteration_frequency(major_iteration_frequency);
  params.set_termination_check_frequency(1);
  params.add_random_projection_seeds(1);
  params.add_random_projection_seeds(2);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  EXPECT_EQ(output.solve_log.iteration_stats().size(), iteration_limit + 1);
  for (const auto& stats : output.solve_log.iteration_stats()) {
    for (const auto& metadata : stats.point_metadata()) {
      // There isn't much we can say about the random projection values, so just
      // check that the right numbers are present.
      EXPECT_THAT(metadata.random_primal_projections(), SizeIs(2));
      EXPECT_THAT(metadata.random_dual_projections(), SizeIs(2));
    }
  }
}

TEST(PrimalDualHybridGradientTest, ProjectInitialPointPrimalBounds) {
  const int iteration_limit = 5;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);

  // The default initial solution (zeros) doesn't satisfy the primal variable
  // bounds. The solver should project it to a valid primal solution.
  SolverResult output =
      PrimalDualHybridGradient(SmallInitializationLp(), params);
  EXPECT_EQ(output.solve_log.iteration_stats().size(), iteration_limit + 1);
  EXPECT_GT(output.primal_solution[0], 0.0);
}

TEST(PrimalDualHybridGradientTest, ProjectInitialPointDualBounds) {
  const int iteration_limit = 5;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  // This initial solution doesn't satisfy the dual variable bounds. The solver
  // should project it to a valid dual solution.
  const PrimalAndDualSolution initial_solution = {
      .primal_solution = Eigen::VectorXd{{1.0, 0.0}},
      .dual_solution = -VectorXd::Ones(2),
  };
  SolverResult output_nonzero_init = PrimalDualHybridGradient(
      SmallInitializationLp(), params, initial_solution);
  EXPECT_EQ(output_nonzero_init.solve_log.iteration_stats().size(),
            iteration_limit + 1);
  EXPECT_LE(output_nonzero_init.dual_solution[0], 0.0);
}

TEST(PrimalDualHybridGradientTest, DetectsProblemWithInconsistentBounds) {
  SolverResult output = PrimalDualHybridGradient(
      SmallInvalidProblemLp(), PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsNonconvexQp) {
  QuadraticProgram qp = TestDiagonalQp1();
  qp.objective_matrix->diagonal()[0] = -1.0;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsProblemWithInconsistentSizes) {
  QuadraticProgram qp = TinyLp();
  qp.objective_vector.resize(0);
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsUncompressedConstraintMatrix) {
  QuadraticProgram qp = TinyLp();
  qp.constraint_matrix.uncompress();
  SolverResult output =
      PrimalDualHybridGradient(std::move(qp), PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsInvalidParameters) {
  PrimalDualHybridGradientParams params;
  params.set_num_threads(0);
  SolverResult output = PrimalDualHybridGradient(TinyLp(), params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PARAMETER);
}

TEST(PrimalDualHybridGradientTest, DetectsNanInConstraintMatrix) {
  QuadraticProgram qp = TestLp();
  qp.constraint_matrix.coeffRef(0, 0) =
      std::numeric_limits<double>::quiet_NaN();
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsExcessiveConstraintMatrix) {
  QuadraticProgram qp = TestLp();
  qp.constraint_matrix.coeffRef(0, 0) = 1.0e60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest,
     DetectsExcessivelySmallColNormConstraintMatrix) {
  QuadraticProgram qp = TestLp();
  qp.constraint_matrix.coeffRef(0, 1) = 1.0e-60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest,
     DetectsExcessivelySmallRowNormConstraintMatrix) {
  QuadraticProgram qp = TestLp();
  qp.constraint_matrix.coeffRef(2, 0) = 1.0e-60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsNanInConstraintBounds) {
  QuadraticProgram qp = TestLp();
  qp.constraint_upper_bounds[1] = std::numeric_limits<double>::quiet_NaN();
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsExcessiveConstraintUpperBound) {
  QuadraticProgram qp = TestLp();
  qp.constraint_upper_bounds[1] = 1.0e60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest,
     ExcessivelySmallConstraintUpperBoundOkWithoutPresolve) {
  QuadraticProgram qp = TestLp();
  qp.constraint_upper_bounds[1] = 1.0e-60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST(PrimalDualHybridGradientTest,
     DetectsExcessivelySmallConstraintUpperBoundWithPresolve) {
  QuadraticProgram qp = TestLp();
  qp.constraint_upper_bounds[1] = 1.0e-60;
  PrimalDualHybridGradientParams params;
  params.mutable_presolve_options()->set_use_glop(true);
  SolverResult output = PrimalDualHybridGradient(qp, params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsExcessiveConstraintLowerBound) {
  QuadraticProgram qp = TestLp();
  qp.constraint_lower_bounds[2] = -1.0e60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest,
     ExcessivelySmallConstraintLowerBoundOkWithoutPresolve) {
  QuadraticProgram qp = TestLp();
  qp.constraint_lower_bounds[2] = 1.0e-60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST(PrimalDualHybridGradientTest,
     DetectsExcessivelySmallConstraintLowerBoundWithPresolve) {
  QuadraticProgram qp = TestLp();
  qp.constraint_lower_bounds[2] = 1.0e-60;
  PrimalDualHybridGradientParams params;
  params.mutable_presolve_options()->set_use_glop(true);
  SolverResult output = PrimalDualHybridGradient(qp, params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsNanInVariableBound) {
  QuadraticProgram qp = TestLp();
  qp.variable_lower_bounds[3] = std::numeric_limits<double>::quiet_NaN();
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsExcessiveVariableBound) {
  QuadraticProgram qp = TestLp();
  qp.variable_lower_bounds[3] = -1.0e60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest,
     ExcessivelySmallVariableBoundOkWithoutPresolve) {
  QuadraticProgram qp = TestLp();
  qp.variable_lower_bounds[1] = 0.0;
  qp.variable_upper_bounds[1] = 1.0e-60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST(PrimalDualHybridGradientTest,
     DetectsExcessivelySmallVariableBoundWithPresolve) {
  QuadraticProgram qp = TestLp();
  qp.variable_lower_bounds[1] = 0.0;
  qp.variable_upper_bounds[1] = 1.0e-60;
  PrimalDualHybridGradientParams params;
  params.mutable_presolve_options()->set_use_glop(true);
  SolverResult output = PrimalDualHybridGradient(qp, params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsNanObjectiveOffset) {
  QuadraticProgram qp = TestLp();
  qp.objective_offset = std::numeric_limits<double>::quiet_NaN();
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsExcessiveObjectiveOffset) {
  QuadraticProgram qp = TestLp();
  qp.objective_offset = -1.0e60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsNanInObjectiveVector) {
  QuadraticProgram qp = TestLp();
  qp.objective_vector[3] = std::numeric_limits<double>::quiet_NaN();
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsExcessiveObjectiveVector) {
  QuadraticProgram qp = TestLp();
  qp.objective_vector[3] = -1.0e60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest,
     ExcessivelySmallObjectiveVectorOkWithoutPresolve) {
  QuadraticProgram qp = TestLp();
  qp.objective_vector[3] = 1.0e-60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST(PrimalDualHybridGradientTest,
     DetectsExcessivelySmallObjectiveVectorWithPresolve) {
  QuadraticProgram qp = TestLp();
  qp.objective_vector[3] = 1.0e-60;
  PrimalDualHybridGradientParams params;
  params.mutable_presolve_options()->set_use_glop(true);
  SolverResult output = PrimalDualHybridGradient(qp, params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsNanInObjectiveMatrix) {
  QuadraticProgram qp = TestDiagonalQp1();
  qp.objective_matrix->diagonal()[0] = std::numeric_limits<double>::quiet_NaN();
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsExcessiveObjectiveMatrix) {
  QuadraticProgram qp = TestDiagonalQp1();
  qp.objective_matrix->diagonal()[0] = 1.0e60;
  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsNanInInitialPrimalSolution) {
  QuadraticProgram qp = TestLp();
  const PrimalAndDualSolution initial_solution = {
      .primal_solution =
          VectorXd{{1.0, std::numeric_limits<double>::quiet_NaN(), 1.0, 1.0}},
      .dual_solution = VectorXd{{1.0, 1.0, 1.0, 1.0}},
  };
  SolverResult output = PrimalDualHybridGradient(
      qp, PrimalDualHybridGradientParams(), initial_solution);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_INITIAL_SOLUTION);
}

TEST(PrimalDualHybridGradientTest,
     DetectsExcessiveValueInInitialPrimalSolution) {
  QuadraticProgram qp = TestLp();
  const PrimalAndDualSolution initial_solution = {
      .primal_solution = VectorXd{{1.0, 1.0e100, 1.0, 1.0}},
      .dual_solution = VectorXd{{1.0, 1.0, 1.0, 1.0}},
  };

  SolverResult output = PrimalDualHybridGradient(
      qp, PrimalDualHybridGradientParams(), initial_solution);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_INITIAL_SOLUTION);
}

TEST(PrimalDualHybridGradientTest, DetectsIncorrectSizeInitialPrimalSolution) {
  QuadraticProgram qp = TestLp();
  const PrimalAndDualSolution initial_solution = {
      .primal_solution = VectorXd{{1.0, 1.0, 1.0}},
      .dual_solution = VectorXd{{1.0, 1.0, 1.0, 1.0}},
  };

  SolverResult output = PrimalDualHybridGradient(
      qp, PrimalDualHybridGradientParams(), initial_solution);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_INITIAL_SOLUTION);
}

TEST(PrimalDualHybridGradientTest, DetectsNanInInitialDualSolution) {
  QuadraticProgram qp = TestLp();
  const PrimalAndDualSolution initial_solution = {
      .primal_solution = VectorXd{{1.0, 1.0, 1.0, 1.0}},
      .dual_solution =
          VectorXd{{1.0, std::numeric_limits<double>::quiet_NaN(), 1.0, 1.0}},
  };
  SolverResult output = PrimalDualHybridGradient(
      qp, PrimalDualHybridGradientParams(), initial_solution);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_INITIAL_SOLUTION);
}

TEST(PrimalDualHybridGradientTest, DetectsExcessiveValueInInitialDualSolution) {
  QuadraticProgram qp = TestLp();
  const PrimalAndDualSolution initial_solution = {
      .primal_solution = VectorXd{{1.0, 1.0, 1.0, 1.0}},
      .dual_solution = VectorXd{{1.0, 1.0e100, 1.0, 1.0}},
  };

  SolverResult output = PrimalDualHybridGradient(
      qp, PrimalDualHybridGradientParams(), initial_solution);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_INITIAL_SOLUTION);
}

TEST(PrimalDualHybridGradientTest, DetectsIncorrectSizeInitialDualSolution) {
  QuadraticProgram qp = TestLp();
  const PrimalAndDualSolution initial_solution = {
      .primal_solution = VectorXd{{1.0, 1.0, 1.0, 1.0}},
      .dual_solution = VectorXd{{1.0, 1.0, 1.0}},
  };

  SolverResult output = PrimalDualHybridGradient(
      qp, PrimalDualHybridGradientParams(), initial_solution);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_INITIAL_SOLUTION);
}

TEST(PrimalDualHybridGradientTest, DetectsZeroObjectiveScalingFactor) {
  QuadraticProgram qp = TestLp();
  qp.objective_scaling_factor = 0.0;

  SolverResult output =
      PrimalDualHybridGradient(qp, PrimalDualHybridGradientParams());
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PrimalDualHybridGradientTest, DetectsFeasibilityPolishingForQp) {
  QuadraticProgram qp = TestDiagonalQp1();
  PrimalDualHybridGradientParams params;
  params.set_use_feasibility_polishing(true);

  SolverResult output = PrimalDualHybridGradient(qp, params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PARAMETER);
}

TEST(PrimalDualHybridGradientTest, ChecksTerminationAtCorrectFrequency) {
  // `termination_check_frequency` is chosen so that it does not divide the
  // `major_iteration_frequency`.
  const int major_iteration_frequency = 5;
  const int termination_check_frequency = 2;
  const int iteration_limit = 16;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_termination_check_frequency(termination_check_frequency);
  params.set_major_iteration_frequency(major_iteration_frequency);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  ASSERT_EQ(output.solve_log.iteration_count(), iteration_limit);

  std::vector<int> termination_checked_at_iteration;

  for (const auto& stats : output.solve_log.iteration_stats()) {
    if (stats.convergence_information_size() > 0) {
      termination_checked_at_iteration.push_back(stats.iteration_number());
    }
  }

  // Termination is checked on the first iteration, the last iteration, every
  // major iteration, and at the `termination_check_frequency` (with a counter
  // reset on every major iteration).
  EXPECT_THAT(termination_checked_at_iteration,
              ElementsAre(0, 2, 4, 5, 7, 9, 10, 12, 14, 15, 16));
}

TEST(PrimalDualHybridGradientTest, CallsCallback) {
  // `termination_check_frequency` is chosen so that it does not divide the
  // `major_iteration_frequency`.
  const int major_iteration_frequency = 5;
  const int termination_check_frequency = 5;
  const int iteration_limit = 16;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(iteration_limit);
  params.set_termination_check_frequency(termination_check_frequency);
  params.set_major_iteration_frequency(major_iteration_frequency);

  int callback_count = 0;
  SolverResult output = PrimalDualHybridGradient(
      TestLp(), params, /*interrupt_solve=*/nullptr,
      /*message_callback=*/nullptr,
      [&callback_count](const IterationCallbackInfo& callback_info) {
        ++callback_count;
      });
  ASSERT_EQ(output.solve_log.iteration_count(), iteration_limit);
  // The callback should be called at every termination check, that is, at
  // iterations 0, 5, 10, 15, and 16, and additionally when returning the final
  // solution.
  EXPECT_EQ(callback_count, 6);
}

// Returns the unique solution of `TinyLp`.
PrimalAndDualSolution TinyLpSolution() {
  return PrimalAndDualSolution{
      .primal_solution = Eigen::VectorXd{{1.0, 0.0, 6.0, 2.0}},
      .dual_solution = Eigen::VectorXd{{0.5, 4.0, 0.0}},
  };
}

TEST(PrimalDualHybridGradientTest, WarmStartedAtOptimum) {
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()
      ->mutable_simple_optimality_criteria()
      ->set_eps_optimal_absolute(1.0e-10);

  SolverResult output =
      PrimalDualHybridGradient(TinyLp(), params, TinyLpSolution());
  // The solver should run a termination check at iteration 0 and determine that
  // the initial point is the solution of the problem.
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({1, 0, 6, 2}, 1.0e-10));
  EXPECT_THAT(output.dual_solution,
              EigenArrayNear<double>({0.5, 4.0, 0.0}, 1.0e-10));
  EXPECT_LE(output.solve_log.iteration_count(), 0);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST(PrimalDualHybridGradientTest, NoMovementWhenStartedAtOptimum) {
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  // Disable scaling because round-off causes a small amount of movement on the
  // first iteration.
  params.set_l2_norm_rescaling(false);
  params.set_l_inf_ruiz_iterations(0);

  const PrimalAndDualSolution initial_solution = TinyLpSolution();
  SolverResult output =
      PrimalDualHybridGradient(TinyLp(), params, initial_solution);
  EXPECT_THAT(output.primal_solution, ElementsAre(1, 0, 6, 2));
  EXPECT_THAT(output.dual_solution, ElementsAre(0.5, 4.0, 0.0));
  EXPECT_LE(output.solve_log.iteration_count(), 1);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST(PrimalDualHybridGradientTest, EmptyQp) {
  MPModelProto proto;
  absl::StatusOr<QuadraticProgram> qp =
      QpFromMpModelProto(proto, /*relax_integer_variables=*/false);
  ASSERT_TRUE(qp.ok()) << qp.status();
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  SolverResult output = PrimalDualHybridGradient(*qp, params);
  EXPECT_THAT(output.primal_solution, ElementsAre());
  EXPECT_THAT(output.dual_solution, ElementsAre());
  EXPECT_EQ(output.solve_log.iteration_count(), 0);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST(PrimalDualHybridGradientTest, RespectsInterrupt) {
  PrimalDualHybridGradientParams params;
  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);

  std::atomic<bool> interrupt_solve = true;
  const SolverResult output =
      PrimalDualHybridGradient(TestLp(), params, &interrupt_solve);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INTERRUPTED_BY_USER);
}

TEST(PrimalDualHybridGradientTest, RespectsInterruptFromCallback) {
  PrimalDualHybridGradientParams params;
  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);

  std::atomic<bool> interrupt_solve = false;
  auto callback = [&](const IterationCallbackInfo& info) {
    if (info.iteration_stats.iteration_number() >= 10) {
      interrupt_solve.store(true);
    }
  };

  const SolverResult output =
      PrimalDualHybridGradient(TestLp(), params, &interrupt_solve,
                               /*message_callback=*/nullptr, callback);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INTERRUPTED_BY_USER);
  EXPECT_GE(output.solve_log.iteration_count(), 10);
}

TEST(PrimalDualHybridGradientTest, IgnoresFalseInterrupt) {
  PrimalDualHybridGradientParams params;
  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);

  std::atomic<bool> interrupt_solve = false;
  const SolverResult output =
      PrimalDualHybridGradient(TestLp(), params, &interrupt_solve);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT);
}

TEST(PrimalDualHybridGradientTest, HugeNumThreads) {
  const int iteration_upperbound = 10;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(
      iteration_upperbound);
  params.set_num_threads(1'000'000'000);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_ITERATION_LIMIT);
}

TEST(PrimalDualHybridGradientTest, HugeNumShards) {
  const int iteration_upperbound = 10;
  PrimalDualHybridGradientParams params = ParamsWithNoLimits();
  params.mutable_termination_criteria()->set_iteration_limit(
      iteration_upperbound);
  params.set_num_shards(1'000'000'000);

  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            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(1.0e-5);
  opt_criteria->set_eps_optimal_primal_residual_relative(0.0);
  opt_criteria->set_eps_optimal_dual_residual_absolute(1.0e-5);
  opt_criteria->set_eps_optimal_dual_residual_relative(0.0);
  opt_criteria->set_eps_optimal_objective_gap_absolute(1.0e-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));
}

// Regression test for b/311455838. Note that this test only fails in debug
// mode, when an infeasible primal variable (from the iterate differences)
// violates presolve's assumptions and triggers a DCHECK() failure.
TEST(PrimalDualHybridGradientTest, IterateDifferenceBoundsInPresolve) {
  // A trivial (but very badly scaled) LP found by fuzzing.
  QuadraticProgram lp(2, 1);
  lp.objective_offset = -3.0e+23;
  lp.objective_vector =
      VectorXd{{2.7369110631344083e-48, -3.0517578125211636e-05}};
  lp.constraint_lower_bounds = VectorXd{{-2.7369110631344083e-48}};
  lp.constraint_upper_bounds = VectorXd{{0}};
  lp.variable_lower_bounds = VectorXd{{1.8446744073709552e+21, -1.0}};
  lp.variable_upper_bounds = VectorXd{{kInfinity, 1.8446744073709552e+21}};
  lp.constraint_matrix.coeffRef(0, 0) = -2.7369110631344083e-48;
  lp.constraint_matrix.coeffRef(0, 1) = 1.0;
  lp.constraint_matrix.makeCompressed();

  PrimalDualHybridGradientParams params;
  params.mutable_termination_criteria()->set_iteration_limit(40);
  auto& presolve_options = *params.mutable_presolve_options();
  presolve_options.set_use_glop(true);
  presolve_options.mutable_glop_parameters()->set_solve_dual_problem(
      operations_research::glop::GlopParameters::LET_SOLVER_DECIDE);
  presolve_options.mutable_glop_parameters()->set_dualizer_threshold(
      1.3574141825331e-312);
  params.set_infinite_constraint_bound_threshold(1.34785525461908e-312);

  SolverResult output = PrimalDualHybridGradient(lp, params);
  EXPECT_THAT(
      output.solve_log.termination_reason(),
      AnyOf(TERMINATION_REASON_ITERATION_LIMIT, TERMINATION_REASON_OPTIMAL));
}

// `FeasibilityPolishingTest` sets `params_` for feasibility polishing, and to
// avoid PDLP features that disrupt a simple analysis of the performance with
// and without feasibility polishing (primal weight adjustments, dynamic step
// sizes, and restarts). It also sets termination criteria with a relaxed
// objective gap tolerance.
class FeasibilityPolishingTest : public ::testing::Test {
 protected:
  FeasibilityPolishingTest() {
    params_.set_linesearch_rule(
        PrimalDualHybridGradientParams::CONSTANT_STEP_SIZE_RULE);
    auto* opt_criteria = params_.mutable_termination_criteria()
                             ->mutable_detailed_optimality_criteria();
    opt_criteria->set_eps_optimal_primal_residual_absolute(1.0e-6);
    opt_criteria->set_eps_optimal_primal_residual_relative(1.0e-6);
    opt_criteria->set_eps_optimal_dual_residual_absolute(1.0e-6);
    opt_criteria->set_eps_optimal_dual_residual_relative(1.0e-6);
    opt_criteria->set_eps_optimal_objective_gap_absolute(1.0e-2);
    opt_criteria->set_eps_optimal_objective_gap_relative(1.0e-2);
    params_.mutable_termination_criteria()->set_iteration_limit(500);
    params_.set_handle_some_primal_gradients_on_finite_bounds_as_residuals(
        false);
    params_.set_primal_weight_update_smoothing(0.0);
    params_.set_use_feasibility_polishing(true);
  }

  PrimalDualHybridGradientParams params_;
};

// `FeasibilityPolishingPrimalTest` contains a very simple `lp_` that requires
// roughly 2500 iterations without feasibility polishing because the dual is
// slow to converge, but with feasibility polishing solves the first time it
// attempts polishing.
class FeasibilityPolishingPrimalTest : public FeasibilityPolishingTest {
 protected:
  FeasibilityPolishingPrimalTest() : lp_(2, 1) {
    // min x_1 + 1.001 x_2 s.t. x_1 + x_2 = 1
    lp_.objective_offset = 0;
    lp_.objective_vector = VectorXd{{1.0, 1.001}};
    lp_.variable_lower_bounds = VectorXd{{-1, -1}};
    lp_.variable_upper_bounds = VectorXd{{kInfinity, kInfinity}};
    lp_.constraint_lower_bounds = VectorXd{{1}};
    lp_.constraint_upper_bounds = VectorXd{{1}};
    lp_.constraint_matrix.coeffRef(0, 0) = 1.0;
    lp_.constraint_matrix.coeffRef(0, 1) = 1.0;
    lp_.constraint_matrix.makeCompressed();
  }

  QuadraticProgram lp_;
};

// `FeasibilityPolishingDualTest` contains the dual of the lp from
// `FeasibilityPolishingPrimalTest` with the lower bounds on the variables
// changed to zero. It requires roughly 2500 iterations without feasibility
// polishing because the primal is slow to converge, but with feasibility
// polishing solves the first time it attempts polishing.
class FeasibilityPolishingDualTest : public FeasibilityPolishingTest {
 protected:
  FeasibilityPolishingDualTest() : lp_(1, 2) {
    // min -y s.t. y <= 1, y <= 1.001, y free
    lp_.objective_offset = 0;
    lp_.objective_vector = VectorXd{{-1.0}};
    lp_.variable_lower_bounds = VectorXd{{-kInfinity}};
    lp_.variable_upper_bounds = VectorXd{{kInfinity}};
    lp_.constraint_lower_bounds = VectorXd{{-kInfinity, -kInfinity}};
    lp_.constraint_upper_bounds = VectorXd{{1, 1.001}};
    lp_.constraint_matrix.coeffRef(0, 0) = 1.0;
    lp_.constraint_matrix.coeffRef(1, 0) = 1.0;
    lp_.constraint_matrix.makeCompressed();
  }

  QuadraticProgram lp_;
};

TEST_F(FeasibilityPolishingPrimalTest, FeasibilityPolishingSolvesFaster) {
  // Without feasibility polishing, PDHG requires roughly 2500 iterations.
  params_.set_use_feasibility_polishing(false);
  SolverResult base_output = PrimalDualHybridGradient(lp_, params_);
  EXPECT_EQ(base_output.solve_log.termination_reason(),
            TERMINATION_REASON_ITERATION_LIMIT);

  // With feasibility polishing, PDHG solves the first time it attempts
  // polishing.
  params_.set_use_feasibility_polishing(true);
  SolverResult output = PrimalDualHybridGradient(lp_, params_);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST_F(FeasibilityPolishingDualTest, FeasibilityPolishingSolvesFaster) {
  // Without feasibility polishing, PDHG requires roughly 2500 iterations.
  params_.set_use_feasibility_polishing(false);
  SolverResult base_output = PrimalDualHybridGradient(lp_, params_);
  EXPECT_EQ(base_output.solve_log.termination_reason(),
            TERMINATION_REASON_ITERATION_LIMIT);

  // With feasibility polishing, PDHG solves the first time it attempts
  // polishing.
  params_.set_use_feasibility_polishing(true);
  SolverResult output = PrimalDualHybridGradient(lp_, params_);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST_F(FeasibilityPolishingPrimalTest, FeasibilityPolishingFindsValidSolution) {
  SolverResult output = PrimalDualHybridGradient(lp_, params_);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
  VerifyObjectiveValues(output, 1.0, 1.0e-2);
  ASSERT_EQ(output.primal_solution.size(), 2);
  // The optimal solution is {1, 0}. However, the relaxed optimality tolerance
  // accepts many solutions, so just check feasibility (in addition to the
  // objective value check above).
  EXPECT_THAT(output.primal_solution[0] + output.primal_solution[1],
              DoubleNear(1.0, 1.0e-6));
  EXPECT_GE(output.primal_solution[0], 0.0);
  EXPECT_GE(output.primal_solution[1], 0.0);

  ASSERT_EQ(output.dual_solution.size(), 1);
  // Dual feasibility.
  EXPECT_LE(output.dual_solution[0], 1.0);
  // Dual optimality.
  EXPECT_GE(output.dual_solution[0], 1.0 - 1e-2);

  EXPECT_THAT(output.reduced_costs,
              EigenArrayNear<double>({1.0 - output.dual_solution[0],
                                      1.001 - output.dual_solution[0]},
                                     1.0e-12));
}

TEST_F(FeasibilityPolishingPrimalTest,
       NoPolishingAfterIterationLimitWhenPolishingAfterLimitsDisabled) {
  // Feasibility polishing would solve the problem the first time it is
  // attempted, which would be at iteration 100.
  params_.set_use_feasibility_polishing(true);
  params_.set_apply_feasibility_polishing_after_limits_reached(false);
  params_.mutable_termination_criteria()->set_iteration_limit(50);
  SolverResult output = PrimalDualHybridGradient(lp_, params_);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_ITERATION_LIMIT);
}

TEST_F(FeasibilityPolishingPrimalTest,
       PolishingAfterIterationLimitWhenPolishingAfterLimitsEnabled) {
  params_.set_use_feasibility_polishing(true);
  params_.set_apply_feasibility_polishing_after_limits_reached(true);
  params_.mutable_termination_criteria()->set_iteration_limit(50);
  SolverResult output = PrimalDualHybridGradient(lp_, params_);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST_F(FeasibilityPolishingPrimalTest,
       PolishingTerminatesAfterIterationLimitWhenPolishingAfterLimitsDisabled) {
  // Feasibility polishing will be triggered at iteration 100. The iteration
  // limit prevents primal polishing from completing.
  params_.set_use_feasibility_polishing(true);
  params_.set_apply_feasibility_polishing_after_limits_reached(false);
  params_.mutable_termination_criteria()->set_iteration_limit(101);
  SolverResult output = PrimalDualHybridGradient(lp_, params_);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_ITERATION_LIMIT);
}

TEST_F(FeasibilityPolishingPrimalTest,
       PolishingContinuesAfterIterationLimitWhenPolishingAfterLimitsEnabled) {
  params_.set_use_feasibility_polishing(true);
  params_.set_apply_feasibility_polishing_after_limits_reached(true);
  params_.mutable_termination_criteria()->set_iteration_limit(101);
  SolverResult output = PrimalDualHybridGradient(lp_, params_);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST_F(FeasibilityPolishingPrimalTest,
       PolishingStopsAfterContinuingAfterIterationLimitWhenNotOptimal) {
  params_.set_use_feasibility_polishing(true);
  params_.set_apply_feasibility_polishing_after_limits_reached(true);
  params_.mutable_termination_criteria()->set_iteration_limit(101);
  auto* opt_criteria = params_.mutable_termination_criteria()
                           ->mutable_detailed_optimality_criteria();
  opt_criteria->set_eps_optimal_primal_residual_absolute(1.0e-16);
  opt_criteria->set_eps_optimal_primal_residual_relative(0.0);
  opt_criteria->set_eps_optimal_dual_residual_absolute(1.0e-16);
  opt_criteria->set_eps_optimal_dual_residual_relative(0.0);
  opt_criteria->set_eps_optimal_objective_gap_absolute(1.0e-16);
  opt_criteria->set_eps_optimal_objective_gap_relative(0.0);
  SolverResult output = PrimalDualHybridGradient(lp_, params_);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_ITERATION_LIMIT);
  // 100 main iterations + at most 12 primal feasibility polishing iterations
  // + at most 12 dual feasibility polishing iterations.
  EXPECT_LE(output.solve_log.iteration_count(), 124);
}

TEST_F(FeasibilityPolishingPrimalTest,
       NoPolishingAfterInterruptWhenPolishingAfterInterruptDisabled) {
  // Feasibility polishing would solve the problem the first time it is
  // attempted, which would be at iteration 100.
  params_.set_use_feasibility_polishing(true);
  params_.set_apply_feasibility_polishing_if_solver_is_interrupted(false);
  std::atomic<bool> interrupt_solve = false;
  auto callback = [&](const IterationCallbackInfo& info) {
    if (info.iteration_stats.iteration_number() >= 50) {
      interrupt_solve.store(true);
    }
  };
  SolverResult output =
      PrimalDualHybridGradient(lp_, params_, &interrupt_solve,
                               /*message_callback=*/nullptr, callback);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INTERRUPTED_BY_USER);
}

TEST_F(FeasibilityPolishingPrimalTest,
       PolishingAfterInterruptWhenPolishingAfterInterruptEnabled) {
  // Feasibility polishing would solve the problem the first time it is
  // attempted, which would be at iteration 100.
  params_.set_use_feasibility_polishing(true);
  params_.set_apply_feasibility_polishing_if_solver_is_interrupted(true);
  std::atomic<bool> interrupt_solve = false;
  auto callback = [&](const IterationCallbackInfo& info) {
    if (info.iteration_stats.iteration_number() >= 50) {
      interrupt_solve.store(true);
    }
  };
  SolverResult output =
      PrimalDualHybridGradient(lp_, params_, &interrupt_solve,
                               /*message_callback=*/nullptr, callback);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST_F(FeasibilityPolishingPrimalTest,
       PolishingTerminatesAfterInterruptWhenPolishingAfterInterruptDisabled) {
  // Feasibility polishing would solve the problem the first time it is
  // attempted, which would be at iteration 100.
  params_.set_use_feasibility_polishing(true);
  params_.set_apply_feasibility_polishing_if_solver_is_interrupted(false);
  std::atomic<bool> interrupt_solve = false;
  auto callback = [&](const IterationCallbackInfo& info) {
    if (info.iteration_type == IterationType::kPrimalFeasibility) {
      interrupt_solve.store(true);
    }
  };
  SolverResult output =
      PrimalDualHybridGradient(lp_, params_, &interrupt_solve,
                               /*message_callback=*/nullptr, callback);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INTERRUPTED_BY_USER);
}

TEST_F(FeasibilityPolishingPrimalTest,
       PolishingContinuesAfterInterruptWhenPolishingAfterInterruptEnabled) {
  // Feasibility polishing would solve the problem the first time it is
  // attempted, which would be at iteration 100.
  params_.set_use_feasibility_polishing(true);
  params_.set_apply_feasibility_polishing_if_solver_is_interrupted(true);
  std::atomic<bool> interrupt_solve = false;
  auto callback = [&](const IterationCallbackInfo& info) {
    if (info.iteration_type == IterationType::kPrimalFeasibility) {
      interrupt_solve.store(true);
    }
  };
  SolverResult output =
      PrimalDualHybridGradient(lp_, params_, &interrupt_solve,
                               /*message_callback=*/nullptr, callback);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
}

TEST_F(FeasibilityPolishingPrimalTest, FeasibilityPolishingDetailsInLog) {
  SolverResult output = PrimalDualHybridGradient(lp_, params_);

  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);

  int primal_optimal_count = 0;
  int dual_optimal_count = 0;
  for (const auto& phase : output.solve_log.feasibility_polishing_details()) {
    switch (phase.polishing_phase_type()) {
      case POLISHING_PHASE_TYPE_PRIMAL_FEASIBILITY: {
        if (phase.termination_reason() == TERMINATION_REASON_OPTIMAL) {
          ++primal_optimal_count;
        }
        break;
      }
      case POLISHING_PHASE_TYPE_DUAL_FEASIBILITY: {
        if (phase.termination_reason() == TERMINATION_REASON_OPTIMAL) {
          ++dual_optimal_count;
        }
        break;
      }
      default: {
        ADD_FAILURE() << "Unexpected polishing_phase_type: "
                      << phase.polishing_phase_type();
        break;
      }
    }
    EXPECT_TRUE(phase.has_solution_stats());
    EXPECT_EQ(phase.solution_stats().iteration_number(),
              phase.iteration_count());
  }
  EXPECT_GE(primal_optimal_count, 1);
  EXPECT_GE(dual_optimal_count, 1);
}

TEST_F(FeasibilityPolishingPrimalTest,
       FeasibilityPolishingUsesInfiniteTolerances) {
  SolverResult output = PrimalDualHybridGradient(lp_, params_);

  for (const auto& phase : output.solve_log.feasibility_polishing_details()) {
    switch (phase.polishing_phase_type()) {
      case POLISHING_PHASE_TYPE_PRIMAL_FEASIBILITY: {
        EXPECT_TRUE(phase.params()
                        .termination_criteria()
                        .has_detailed_optimality_criteria());
        const auto& criteria = phase.params()
                                   .termination_criteria()
                                   .detailed_optimality_criteria();
        EXPECT_EQ(criteria.eps_optimal_objective_gap_absolute(), kInfinity);
        EXPECT_EQ(criteria.eps_optimal_objective_gap_relative(), kInfinity);
        EXPECT_EQ(criteria.eps_optimal_dual_residual_absolute(), kInfinity);
        EXPECT_EQ(criteria.eps_optimal_dual_residual_relative(), kInfinity);
        break;
      }
      case POLISHING_PHASE_TYPE_DUAL_FEASIBILITY: {
        EXPECT_TRUE(phase.params()
                        .termination_criteria()
                        .has_detailed_optimality_criteria());
        const auto& criteria = phase.params()
                                   .termination_criteria()
                                   .detailed_optimality_criteria();
        EXPECT_EQ(criteria.eps_optimal_objective_gap_absolute(), kInfinity);
        EXPECT_EQ(criteria.eps_optimal_objective_gap_relative(), kInfinity);
        EXPECT_EQ(criteria.eps_optimal_primal_residual_absolute(), kInfinity);
        EXPECT_EQ(criteria.eps_optimal_primal_residual_relative(), kInfinity);
        break;
      }
      default: {
        ADD_FAILURE() << "Unexpected polishing_phase_type: "
                      << phase.polishing_phase_type();
        break;
      }
    }
  }
}

TEST_F(FeasibilityPolishingPrimalTest,
       FeasibilityPolishingWorkStatsAreCorrect) {
  params_.set_record_iteration_stats(true);
  SolverResult output = PrimalDualHybridGradient(lp_, params_);

  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);

  int total_feasibility_iterations = 0;
  double total_feasibility_time = 0.0;
  for (const auto& phase : output.solve_log.feasibility_polishing_details()) {
    EXPECT_TRUE(phase.has_solution_stats());
    EXPECT_EQ(phase.solution_stats().iteration_number(),
              phase.iteration_count());
    // For at least one of these, `phase.iteration_count() > 0`, because
    // `total_feasibility_iterations` is verified to be >= 1 below.
    EXPECT_GE(phase.iteration_count(), 0);
    EXPECT_LE(phase.iteration_count(), phase.main_iteration_count() / 4);
    total_feasibility_iterations += phase.iteration_count();
    total_feasibility_time += phase.solve_time_sec();
  }
  ASSERT_THAT(output.solve_log.iteration_stats(), Not(IsEmpty()));
  const auto& last_stats = output.solve_log.iteration_stats(
      output.solve_log.iteration_stats_size() - 1);
  EXPECT_GE(total_feasibility_iterations, 1);
  EXPECT_GE(last_stats.iteration_number(), 1);
  // This checks that `iteration_count()` includes both the main and feasibility
  // iterations, and that `iteration_stats.iteration_number()` does not include
  // work from feasibility iterations.
  EXPECT_EQ(last_stats.iteration_number() + total_feasibility_iterations,
            output.solve_log.iteration_count());
  EXPECT_GT(total_feasibility_time, 0.0);
  EXPECT_GT(last_stats.cumulative_time_sec(), 0.0);
  // This is an approximate check that `solve_time_sec()` includes both the main
  // and feasibility iterations, and that
  // `iteration_stats.cumulative_time_sec()` does not include time from
  // feasibility iterations. We don't expect equality because some clock time
  // can pass while switching between phases.
  EXPECT_LE(total_feasibility_time + last_stats.cumulative_time_sec(),
            output.solve_log.solve_time_sec());
}

TEST_F(FeasibilityPolishingPrimalTest, FeasibilityObjectivesAreZero) {
  params_.set_record_iteration_stats(true);
  SolverResult output = PrimalDualHybridGradient(lp_, params_);

  int primal_convergence_info_count = 0;
  int dual_convergence_info_count = 0;
  for (const auto& phase : output.solve_log.feasibility_polishing_details()) {
    for (const auto& iter_stats : phase.iteration_stats()) {
      for (const auto& convergence : iter_stats.convergence_information()) {
        switch (phase.polishing_phase_type()) {
          case POLISHING_PHASE_TYPE_PRIMAL_FEASIBILITY: {
            ++primal_convergence_info_count;
            EXPECT_EQ(convergence.primal_objective(), 0.0);
            break;
          }
          case POLISHING_PHASE_TYPE_DUAL_FEASIBILITY: {
            ++dual_convergence_info_count;
            EXPECT_EQ(convergence.dual_objective(), 0.0);
            break;
          }
          default: {
            ADD_FAILURE() << "Unexpected polishing_phase_type: "
                          << phase.polishing_phase_type();
            break;
          }
        }
      }
    }
  }

  // These checks both verify that the previous objective tests were actually
  // applied, and that iteration stats (and their convergence info) are logged.
  EXPECT_GE(primal_convergence_info_count, 1);
  EXPECT_GE(dual_convergence_info_count, 1);
}

TEST_F(FeasibilityPolishingPrimalTest, CallsCallbackForAllThreePhases) {
  absl::flat_hash_map<IterationType, int> callback_count;
  bool polishing_termination_is_last = false;
  SolverResult output = PrimalDualHybridGradient(
      lp_, params_, /*interrupt_solve=*/nullptr,
      /*message_callback=*/nullptr,
      [&](const IterationCallbackInfo& callback_info) {
        ++callback_count[callback_info.iteration_type];
        polishing_termination_is_last =
            callback_info.iteration_type ==
            IterationType::kFeasibilityPolishingTermination;
      });
  EXPECT_TRUE(polishing_termination_is_last);
  EXPECT_THAT(
      callback_count,
      UnorderedElementsAre(
          Pair(IterationType::kPrimalFeasibility, Ge(1)),
          Pair(IterationType::kDualFeasibility, Ge(1)),
          Pair(IterationType::kNormal, Ge(1)),
          Pair(IterationType::kFeasibilityPolishingTermination, Ge(1))));
}

TEST_F(FeasibilityPolishingTest, IterationLimitIncludesFeasibilityPhases) {
  params_.mutable_termination_criteria()->set_iteration_limit(300);
  params_.set_record_iteration_stats(true);
  SolverResult output = PrimalDualHybridGradient(TinyLp(), params_);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
  int feasibility_polishing_iterations = 0;
  for (const auto& phase : output.solve_log.feasibility_polishing_details()) {
    feasibility_polishing_iterations += phase.iteration_count();
  }
  // `feasibility_polishing_iterations` should at least include 12 iterations
  // from the first failed primal feasibility polishing attempt.
  EXPECT_GE(feasibility_polishing_iterations, 12);
  ASSERT_THAT(output.solve_log.iteration_stats(), Not(IsEmpty()));
  const IterationStats& last_stats = output.solve_log.iteration_stats(
      output.solve_log.iteration_stats_size() - 1);
  EXPECT_EQ(last_stats.iteration_number() + feasibility_polishing_iterations,
            output.solve_log.iteration_count());
}

// This test doesn't attempt to check what is logged. Rather, it just verifies
// that the code succeeds at each verbosity level. Other than having the
// verbosity level as a parameter, and fixing the other parameters, it is the
// same as `PrimalDualHybridGradientLPTest.Tiny`.
TEST_P(PrimalDualHybridGradientVerbosityTest, TinyLp) {
  const int iteration_upperbound = 300;
  const int verbosity_level = GetParam();
  PrimalDualHybridGradientParams params =
      CreateSolverParams(iteration_upperbound,
                         /*eps_optimal_absolute=*/1.0e-5,
                         /*enable_scaling=*/false,
                         /*num_threads=*/1,
                         /*use_iteration_limit=*/false,
                         /*use_malitsky_pock_linesearch=*/false,
                         /*use_diagonal_qp_trust_region_solver=*/false);
  params.set_major_iteration_frequency(60);
  params.set_verbosity_level(verbosity_level);

  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));
  EXPECT_EQ(output.solve_log.original_problem_stats().num_variables(), 4);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_variables(), 4);
  EXPECT_EQ(output.solve_log.original_problem_stats().num_constraints(), 3);
  EXPECT_LE(output.solve_log.preprocessed_problem_stats().num_constraints(), 3);
}

TEST(PresolveTest, DetectsProblemWithInconsistentBounds) {
  PrimalDualHybridGradientParams params;
  params.mutable_presolve_options()->set_use_glop(true);
  SolverResult output =
      PrimalDualHybridGradient(SmallInvalidProblemLp(), params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INVALID_PROBLEM);
}

TEST(PresolveTest, PresolveSolvesToOptimality) {
  PrimalDualHybridGradientParams params;
  params.mutable_presolve_options()->set_use_glop(true);
  SolverResult output = PrimalDualHybridGradient(TestLp(), params);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
  EXPECT_EQ(output.solve_log.iteration_count(), 0);
  ASSERT_EQ(output.solve_log.solution_type(), POINT_TYPE_PRESOLVER_SOLUTION);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({-1, 8, 1, 2.5}, 1.0e-10));
  EXPECT_THAT(output.dual_solution,
              EigenArrayNear<double>({-2, 0, 2.375, 2.0 / 3}, 1.0e-10));
  const auto convergence_info = GetConvergenceInformation(
      output.solve_log.solution_stats(), POINT_TYPE_PRESOLVER_SOLUTION);
  ASSERT_TRUE(convergence_info.has_value());
  EXPECT_EQ(convergence_info->candidate_type(), POINT_TYPE_PRESOLVER_SOLUTION);
  EXPECT_DOUBLE_EQ(convergence_info->primal_objective(), -34.0);
  EXPECT_DOUBLE_EQ(convergence_info->dual_objective(), -34.0);
  EXPECT_DOUBLE_EQ(convergence_info->corrected_dual_objective(), -34.0);
  EXPECT_DOUBLE_EQ(convergence_info->l_inf_primal_residual(), 0.0);
  EXPECT_DOUBLE_EQ(convergence_info->l2_primal_residual(), 0.0);
  EXPECT_DOUBLE_EQ(convergence_info->l_inf_dual_residual(), 0.0);
  EXPECT_DOUBLE_EQ(convergence_info->l2_dual_residual(), 0.0);
  EXPECT_DOUBLE_EQ(convergence_info->l_inf_primal_variable(), 8.0);
  EXPECT_DOUBLE_EQ(convergence_info->l2_primal_variable(), 8.5);
  EXPECT_DOUBLE_EQ(convergence_info->l_inf_dual_variable(), 2.375);
  EXPECT_THAT(convergence_info->l2_dual_variable(),
              DoubleNear(3.17569983538, 1.0e-10));
}

TEST(PresolveTest, SolvesWithGlopScalingEnabled) {
  PrimalDualHybridGradientParams params;
  params.mutable_presolve_options()->set_use_glop(true);
  params.mutable_presolve_options()->mutable_glop_parameters()->set_use_scaling(
      true);
  params.mutable_presolve_options()
      ->mutable_glop_parameters()
      ->set_use_preprocessing(false);
  QuadraticProgram test_lp = TestLp();
  // Change `objective_scaling_factor`, so that glop's presolve scaling will
  // react.
  test_lp.objective_scaling_factor = 0.1;
  SolverResult output = PrimalDualHybridGradient(test_lp, params);
  VerifyTerminationReasonAndIterationCount(params, output,
                                           /*use_iteration_limit=*/false);
  VerifyObjectiveValues(output, -3.4, 1.0e-6);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({-1, 8, 1, 2.5}, 1.0e-4));
  EXPECT_THAT(output.dual_solution,
              EigenArrayNear<double>({-2, 0, 2.375, 2.0 / 3}, 1.0e-4));
}

TEST(PresolveTest, PresolveInfeasible) {
  PrimalDualHybridGradientParams params;
  params.mutable_presolve_options()->set_use_glop(true);
  SolverResult output =
      PrimalDualHybridGradient(SmallPrimalInfeasibleLp(), params);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_PRIMAL_OR_DUAL_INFEASIBLE);
  EXPECT_EQ(output.solve_log.solution_type(), POINT_TYPE_PRESOLVER_SOLUTION);
  EXPECT_EQ(output.solve_log.iteration_count(), 0);
}

TEST(PresolveTest, WarnsInitialSolution) {
  PrimalAndDualSolution initial_solution;
  initial_solution.primal_solution.setZero(6);
  initial_solution.dual_solution.setZero(3);
  PrimalDualHybridGradientParams params;
  params.mutable_presolve_options()->set_use_glop(true);
  SolverResult output = PrimalDualHybridGradient(CorrelationClusteringStarLp(),
                                                 params, initial_solution);
}

TEST(PresolveTest, WarnsInitialSolutionViaCallback) {
  PrimalAndDualSolution initial_solution;
  initial_solution.primal_solution.setZero(6);
  initial_solution.dual_solution.setZero(3);
  PrimalDualHybridGradientParams params;
  params.mutable_presolve_options()->set_use_glop(true);
  std::vector<std::string> message_output;
  SolverResult output = PrimalDualHybridGradient(
      CorrelationClusteringStarLp(), params, initial_solution,
      /*interrupt_solve=*/nullptr,
      [&message_output](const std::string& message) {
        message_output.push_back(message);
      });
  EXPECT_THAT(message_output, Contains(AllOf(HasSubstr("initial solution"),
                                             HasSubstr("presolve"))));
}

TEST_P(PresolveDualScalingTest, Dualize) {
  auto [dualize, negate_and_scale_objective] = GetParam();
  PrimalDualHybridGradientParams params;
  params.mutable_presolve_options()->set_use_glop(true);
  params.mutable_presolve_options()
      ->mutable_glop_parameters()
      ->set_solve_dual_problem(dualize ? glop::GlopParameters::ALWAYS_DO
                                       : glop::GlopParameters::NEVER_DO);
  params.mutable_termination_criteria()->set_iteration_limit(1000);
  QuadraticProgram qp = CorrelationClusteringStarLp();
  if (negate_and_scale_objective) {
    qp.objective_scaling_factor = -3;
  }
  SolverResult output = PrimalDualHybridGradient(qp, params);
  EXPECT_EQ(output.solve_log.termination_reason(), TERMINATION_REASON_OPTIMAL);
  EXPECT_GT(output.solve_log.iteration_count(), 0);
  EXPECT_THAT(output.primal_solution,
              EigenArrayNear<double>({0.5, 0.5, 0.5, 0.0, 0.0, 0.0}, 1.0e-10));
  EXPECT_THAT(output.dual_solution,
              EigenArrayNear<double>({0.5, 0.5, 0.5}, 1.0e-10));
}

TEST(PresolveTest, PresolveParametersAreUsed) {
  QuadraticProgram qp = CorrelationClusteringStarLp();
  PrimalDualHybridGradientParams params;
  params.mutable_termination_criteria()->set_iteration_limit(0);
  params.mutable_presolve_options()->set_use_glop(true);
  SolverResult output_1 = PrimalDualHybridGradient(qp, params);
  params.mutable_presolve_options()
      ->mutable_glop_parameters()
      ->set_solve_dual_problem(glop::GlopParameters::ALWAYS_DO);
  SolverResult output_2 = PrimalDualHybridGradient(qp, params);
  EXPECT_NE(output_1.solve_log.preprocessed_problem_stats().num_variables(),
            output_2.solve_log.preprocessed_problem_stats().num_variables());
  EXPECT_NE(output_1.solve_log.preprocessed_problem_stats().num_constraints(),
            output_2.solve_log.preprocessed_problem_stats().num_constraints());
}

TEST(ComputeStatusesTest, AtOptimum) {
  QuadraticProgram lp = TestLp();
  const PrimalAndDualSolution solution = {
      .primal_solution = Eigen::VectorXd{{-1, 8, 1, 2.5}},
      .dual_solution = Eigen::VectorXd{{-2, 0, 2.375, 2.0 / 3}},
  };
  glop::ProblemSolution glop_solution = internal::ComputeStatuses(lp, solution);
  EXPECT_THAT(
      glop_solution.constraint_statuses,
      ElementsAre(ConstraintStatus::FIXED_VALUE, ConstraintStatus::BASIC,
                  ConstraintStatus::AT_LOWER_BOUND,
                  ConstraintStatus::AT_LOWER_BOUND));
  EXPECT_THAT(
      glop_solution.variable_statuses,
      ElementsAre(VariableStatus::BASIC, VariableStatus::BASIC,
                  VariableStatus::BASIC, VariableStatus::AT_LOWER_BOUND));
}

TEST(ComputeStatusesTest, CoverMoreCases) {
  QuadraticProgram lp = TestLp();
  lp.variable_upper_bounds[3] = 2.5;
  lp.variable_upper_bounds[1] = 8;
  const PrimalAndDualSolution solution = {
      .primal_solution = Eigen::VectorXd{{-1, 8, 1, 2.5}},
      .dual_solution = Eigen::VectorXd{{-2, 0, 2.375, -1}},
  };
  glop::ProblemSolution glop_solution = internal::ComputeStatuses(lp, solution);
  EXPECT_THAT(
      glop_solution.constraint_statuses,
      ElementsAre(ConstraintStatus::FIXED_VALUE, ConstraintStatus::BASIC,
                  ConstraintStatus::AT_LOWER_BOUND,
                  ConstraintStatus::AT_UPPER_BOUND));
  EXPECT_THAT(glop_solution.variable_statuses,
              ElementsAre(VariableStatus::BASIC, VariableStatus::AT_UPPER_BOUND,
                          VariableStatus::BASIC, VariableStatus::FIXED_VALUE));
}

}  // namespace
}  // namespace operations_research::pdlp