primal_dual_hybrid_gradient_test.cc

Go to the documentation of this file.

// Copyright 2010-2022 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//     http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
 
#include "ortools/pdlp/primal_dual_hybrid_gradient.h"
 
#include <algorithm>
#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/log/check.h"
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "absl/strings/str_cat.h"
#include "gmock/gmock.h"
#include "gtest/gtest.h"
#include "ortools/base/logging.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 ::operations_research::glop::ConstraintStatus;
using ::operations_research::glop::VariableStatus;
using ::testing::_;
using ::testing::AnyNumber;
using ::testing::AnyOf;
using ::testing::DoubleNear;
using ::testing::ElementsAre;
using ::testing::Eq;
using ::testing::HasSubstr;
using ::testing::IsEmpty;
using ::testing::Not;
using ::testing::SizeIs;
 
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);
  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>({std::max(dual[1] - dual[0], 0.0),
                                      std::max(dual[0] - dual[1], 0.0)},
                                     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, 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);
 
  PrimalAndDualSolution initial_solution;
  initial_solution.primal_solution.resize(2);
  initial_solution.primal_solution << 0.999, 0.001;
  initial_solution.dual_solution.resize(1);
  initial_solution.dual_solution << -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 << -1, -kInfinity, -2;
  qp.variable_upper_bounds << kInfinity, 4, 10;
  qp.objective_vector << 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 << -1, -kInfinity, -2;
  qp.constraint_upper_bounds << 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 << 1;
  qp.variable_upper_bounds << 1;
  qp.objective_vector << 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 << -1;
  qp.constraint_upper_bounds << -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 number 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.
  PrimalAndDualSolution initial_solution;
  initial_solution.primal_solution = Eigen::VectorXd(2);
  initial_solution.primal_solution << 1.0, 0.0;
  initial_solution.dual_solution = -Eigen::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, 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, 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, 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, DetectsExcessiveVariableBoundGap) {
  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, 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, 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();
  PrimalAndDualSolution initial_solution = {
      .primal_solution =
          Eigen::VectorXd{
              {1.0, std::numeric_limits<double>::quiet_NaN(), 1.0, 1.0}},
      .dual_solution = Eigen::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();
  PrimalAndDualSolution initial_solution = {
      .primal_solution = Eigen::VectorXd{{1.0, 1.0e100, 1.0, 1.0}},
      .dual_solution = Eigen::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();
  PrimalAndDualSolution initial_solution = {
      .primal_solution = Eigen::VectorXd{{1.0, 1.0, 1.0}},
      .dual_solution = Eigen::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();
  PrimalAndDualSolution initial_solution = {
      .primal_solution = Eigen::VectorXd{{1.0, 1.0, 1.0, 1.0}},
      .dual_solution =
          Eigen::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();
  PrimalAndDualSolution initial_solution = {
      .primal_solution = Eigen::VectorXd{{1.0, 1.0, 1.0, 1.0}},
      .dual_solution = Eigen::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();
  PrimalAndDualSolution initial_solution = {
      .primal_solution = Eigen::VectorXd{{1.0, 1.0, 1.0, 1.0}},
      .dual_solution = Eigen::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, 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,
      [&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.
  CHECK_EQ(callback_count, 6);
}
 
// Returns the unique solution of `TinyLp`.
PrimalAndDualSolution TinyLpSolution() {
  PrimalAndDualSolution solution;
  solution.primal_solution.resize(4);
  solution.primal_solution << 1.0, 0.0, 6.0, 2.0;
  solution.dual_solution.resize(3);
  solution.dual_solution << 0.5, 4.0, 0.0;
  return solution;
}
 
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);
 
  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) {
  std::atomic<bool> interrupt_solve;
  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);
 
  interrupt_solve.store(true);
  const SolverResult output =
      PrimalDualHybridGradient(TestLp(), params, &interrupt_solve);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INTERRUPTED_BY_USER);
}
 
TEST(PrimalDualHybridGradientTest, RespectsInterruptFromCallback) {
  std::atomic<bool> interrupt_solve;
  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);
 
  interrupt_solve.store(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, callback);
  EXPECT_EQ(output.solve_log.termination_reason(),
            TERMINATION_REASON_INTERRUPTED_BY_USER);
  EXPECT_GE(output.solve_log.iteration_count(), 10);
}
 
TEST(PrimalDualHybridGradientTest, IgnoresFalseInterrupt) {
  std::atomic<bool> interrupt_solve;
  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);
 
  interrupt_solve.store(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));
}
 
// Verifies that the primal and dual solution satisfy the bounds constraints.
// This function uses ASSERTS rather than EXPECTs because it's used with large
// QPs that would spam the logs otherwise.
void VerifyBoundConstraints(const QuadraticProgram& qp,
                            const Eigen::VectorXd primal_solution,
                            const Eigen::VectorXd dual_solution) {
  for (int64_t i = 0; i < primal_solution.size(); ++i) {
    ASSERT_TRUE(std::isfinite(primal_solution[i]))
        << i << " " << primal_solution[i];
    ASSERT_GE(primal_solution[i], qp.variable_lower_bounds[i]);
    ASSERT_LE(primal_solution[i], qp.variable_upper_bounds[i]);
  }
  for (int i = 0; i < dual_solution.size(); ++i) {
    ASSERT_TRUE(std::isfinite(dual_solution[i]))
        << i << " " << dual_solution[i];
    if (!std::isfinite(qp.constraint_lower_bounds[i])) {
      ASSERT_LE(dual_solution[i], 0);
    }
    if (!std::isfinite(qp.constraint_upper_bounds[i])) {
      ASSERT_GE(dual_solution[i], 0);
    }
  }
}
 
// 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_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();
  PrimalAndDualSolution solution;
  solution.primal_solution.resize(4);
  solution.primal_solution << -1, 8, 1, 2.5;
  solution.dual_solution.resize(4);
  solution.dual_solution << -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;
  PrimalAndDualSolution solution;
  solution.primal_solution.resize(4);
  solution.primal_solution << -1, 8, 1, 2.5;
  solution.dual_solution.resize(4);
  solution.dual_solution << -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