// 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/iteration_stats.h"

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <limits>
#include <optional>
#include <random>
#include <utility>
#include <vector>

#include "Eigen/Core"
#include "Eigen/SparseCore"
#include "absl/log/check.h"
#include "absl/random/distributions.h"
#include "ortools/base/mathutil.h"
#include "ortools/pdlp/quadratic_program.h"
#include "ortools/pdlp/sharded_quadratic_program.h"
#include "ortools/pdlp/sharder.h"
#include "ortools/pdlp/solve_log.pb.h"
#include "ortools/pdlp/solvers.pb.h"

namespace operations_research::pdlp {
namespace {

using ::Eigen::VectorXd;

// `ResidualNorms` contains measures of the infeasibility of a primal or dual
// solution. `objective_correction` is the (additive) adjustment to the
// objective function from the reduced costs. `objective_full_correction` is the
// (additive) adjustment to the objective function if all dual residuals were
// set to zero, while `l_inf_residual`, `l_2_residual`, and
// `l_inf_componentwise_residual` are the L_infinity, L_2, and L_infinity
// (componentwise) norms of the residuals (portions of the primal gradient not
// included in the reduced costs).
struct ResidualNorms {
  double objective_correction;
  double objective_full_correction;
  double l_inf_residual;
  double l_2_residual;
  double l_inf_componentwise_residual;
};

// Computes norms of the primal residual infeasibilities (b - A x) of the
// unscaled problem. Note the primal residuals of the unscaled problem are equal
// to those of the scaled problem divided by `row_scaling_vec`. Does not perform
// any corrections (so the returned `.objective_correction == 0` and
// `.objective_full_correction == 0`). `sharded_qp` is assumed to be the scaled
// problem. If `use_homogeneous_constraint_bounds` is set to true the residuals
// are computed with all finite bounds mapped to zero.
// NOTE: `componentwise_residual_offset` only affects the value of
// `l_inf_componentwise_residual` in the returned `ResidualNorms`.
ResidualNorms PrimalResidualNorms(
    const ShardedQuadraticProgram& sharded_qp, const VectorXd& row_scaling_vec,
    const VectorXd& scaled_primal_solution,
    const double componentwise_residual_offset,
    bool use_homogeneous_constraint_bounds = false) {
  const QuadraticProgram& qp = sharded_qp.Qp();
  CHECK_EQ(row_scaling_vec.size(), sharded_qp.DualSize());
  CHECK_EQ(scaled_primal_solution.size(), sharded_qp.PrimalSize());

  VectorXd primal_product = TransposedMatrixVectorProduct(
      sharded_qp.TransposedConstraintMatrix(), scaled_primal_solution,
      sharded_qp.TransposedConstraintMatrixSharder());
  VectorXd local_l_inf_residual(sharded_qp.DualSharder().NumShards());
  VectorXd local_sumsq_residual(sharded_qp.DualSharder().NumShards());
  VectorXd local_l_inf_componentwise_residual(
      sharded_qp.DualSharder().NumShards());
  sharded_qp.DualSharder().ParallelForEachShard(
      [&](const Sharder::Shard& shard) {
        const auto lower_bound_shard = shard(qp.constraint_lower_bounds);
        const auto upper_bound_shard = shard(qp.constraint_upper_bounds);
        const auto row_scaling_shard = shard(row_scaling_vec);
        const auto primal_product_shard = shard(primal_product);
        double l_inf_residual = 0.0;
        double sumsq_residual = 0.0;
        double l_inf_componentwise_residual = 0.0;
        for (int64_t i = 0; i < primal_product_shard.size(); ++i) {
          const double upper_bound = (use_homogeneous_constraint_bounds &&
                                      std::isfinite(upper_bound_shard[i]))
                                         ? 0.0
                                         : upper_bound_shard[i];
          const double lower_bound = (use_homogeneous_constraint_bounds &&
                                      std::isfinite(lower_bound_shard[i]))
                                         ? 0.0
                                         : lower_bound_shard[i];
          double scaled_residual = 0.0;
          double residual_bound = 0.0;
          if (primal_product_shard[i] > upper_bound) {
            scaled_residual = primal_product_shard[i] - upper_bound;
            residual_bound = upper_bound;
          } else if (primal_product_shard[i] < lower_bound) {
            scaled_residual = lower_bound - primal_product_shard[i];
            residual_bound = lower_bound;
          }
          const double residual = scaled_residual / row_scaling_shard[i];
          l_inf_residual = std::max(l_inf_residual, residual);
          sumsq_residual += residual * residual;
          // Special case: ignore `residual` if == 0, to avoid NaN if offset and
          // bound are both zero.
          if (residual > 0.0) {
            l_inf_componentwise_residual = std::max(
                l_inf_componentwise_residual,
                residual / (componentwise_residual_offset +
                            std::abs(residual_bound / row_scaling_shard[i])));
          }
        }
        local_l_inf_residual[shard.Index()] = l_inf_residual;
        local_sumsq_residual[shard.Index()] = sumsq_residual;
        local_l_inf_componentwise_residual[shard.Index()] =
            l_inf_componentwise_residual;
      });
  return ResidualNorms{
      .objective_correction = 0.0,
      .objective_full_correction = 0.0,
      .l_inf_residual = local_l_inf_residual.lpNorm<Eigen::Infinity>(),
      .l_2_residual = std::sqrt(local_sumsq_residual.sum()),
      .l_inf_componentwise_residual =
          local_l_inf_componentwise_residual.lpNorm<Eigen::Infinity>(),
  };
}

bool TreatVariableBoundAsFinite(const PrimalDualHybridGradientParams& params,
                                double primal_value, double bound) {
  if (params.handle_some_primal_gradients_on_finite_bounds_as_residuals()) {
    // Note that this test is always false if `bound` is infinite.
    return std::abs(primal_value - bound) <= std::abs(primal_value);
  } else {
    return std::isfinite(bound);
  }
}

struct VariableBounds {
  double lower_bound;
  double upper_bound;
};

VariableBounds EffectiveVariableBounds(
    const PrimalDualHybridGradientParams& params, double primal_value,
    double lower_bound, double upper_bound) {
  return {.lower_bound =
              TreatVariableBoundAsFinite(params, primal_value, lower_bound)
                  ? lower_bound
                  : -std::numeric_limits<double>::infinity(),
          .upper_bound =
              TreatVariableBoundAsFinite(params, primal_value, upper_bound)
                  ? upper_bound
                  : std::numeric_limits<double>::infinity()};
}

double VariableBoundForDualObjective(double primal_gradient,
                                     const VariableBounds& bounds) {
  const double primary_bound =
      primal_gradient >= 0.0 ? bounds.lower_bound : bounds.upper_bound;
  const double secondary_bound =
      primal_gradient >= 0.0 ? bounds.upper_bound : bounds.lower_bound;
  if (std::isfinite(primary_bound)) {
    return primary_bound;
  } else if (std::isfinite(secondary_bound)) {
    return secondary_bound;
  } else {
    return 0.0;
  }
}

// Computes norms of the dual residuals and reduced costs of the unscaled
// problem. Note the primal gradient of the unscaled problem is equal to
// `scaled_primal_gradient` divided by `col_scaling_vec`. `sharded_qp` is
// assumed to be the scaled problem. See
// https://developers.google.com/optimization/lp/pdlp_math and the documentation
// for `PrimalDualHybridGradientParams::
// handle_some_primal_gradients_on_finite_bounds_as_residuals` for details and
// notation.
// NOTE: `componentwise_residual_offset` only affects the value of
// `l_inf_componentwise_residual` in the returned `ResidualNorms`.
ResidualNorms DualResidualNorms(const PrimalDualHybridGradientParams& params,
                                const ShardedQuadraticProgram& sharded_qp,
                                const VectorXd& col_scaling_vec,
                                const VectorXd& scaled_primal_solution,
                                const VectorXd& scaled_primal_gradient,
                                const double componentwise_residual_offset) {
  const QuadraticProgram& qp = sharded_qp.Qp();
  CHECK_EQ(col_scaling_vec.size(), sharded_qp.PrimalSize());
  CHECK_EQ(scaled_primal_gradient.size(), sharded_qp.PrimalSize());
  VectorXd local_dual_correction(sharded_qp.PrimalSharder().NumShards());
  VectorXd local_dual_full_correction(sharded_qp.PrimalSharder().NumShards());
  VectorXd local_l_inf_residual(sharded_qp.PrimalSharder().NumShards());
  VectorXd local_sumsq_residual(sharded_qp.PrimalSharder().NumShards());
  VectorXd local_l_inf_componentwise_residual(
      sharded_qp.PrimalSharder().NumShards());
  sharded_qp.PrimalSharder().ParallelForEachShard(
      [&](const Sharder::Shard& shard) {
        const auto lower_bound_shard = shard(qp.variable_lower_bounds);
        const auto upper_bound_shard = shard(qp.variable_upper_bounds);
        const auto primal_gradient_shard = shard(scaled_primal_gradient);
        const auto col_scaling_shard = shard(col_scaling_vec);
        const auto primal_solution_shard = shard(scaled_primal_solution);
        const auto objective_shard = shard(qp.objective_vector);
        double dual_correction = 0.0;
        double dual_full_correction = 0.0;
        double l_inf_residual = 0.0;
        double sumsq_residual = 0.0;
        double l_inf_componentwise_residual = 0.0;
        for (int64_t i = 0; i < primal_gradient_shard.size(); ++i) {
          // The corrections use the scaled values because
          // unscaled_lower_bound = lower_bound * scale and
          // unscaled_primal_gradient = primal_gradient / scale, so the scales
          // cancel out.
          if (primal_gradient_shard[i] == 0.0) continue;
          const double upper_bound = upper_bound_shard[i];
          const double lower_bound = lower_bound_shard[i];
          const double bound_for_rc =
              primal_gradient_shard[i] > 0.0 ? lower_bound : upper_bound;
          dual_full_correction += bound_for_rc * primal_gradient_shard[i];
          VariableBounds effective_bounds = EffectiveVariableBounds(
              params, primal_solution_shard[i], lower_bound, upper_bound);
          // The dual correction (using the appropriate bound) is applied even
          // if the gradient is handled as a residual, so that the dual
          // objective is convex.
          dual_correction += VariableBoundForDualObjective(
                                 primal_gradient_shard[i], effective_bounds) *
                             primal_gradient_shard[i];
          const double effective_bound_for_residual =
              primal_gradient_shard[i] > 0.0 ? effective_bounds.lower_bound
                                             : effective_bounds.upper_bound;
          if (std::isinf(effective_bound_for_residual)) {
            const double scaled_residual = std::abs(primal_gradient_shard[i]);
            const double residual = scaled_residual / col_scaling_shard[i];
            l_inf_residual = std::max(l_inf_residual, residual);
            sumsq_residual += residual * residual;
            // Special case: ignore `residual` if == 0, to avoid NaN if offset
            // and objective are both zero.
            if (residual > 0.0) {
              l_inf_componentwise_residual = std::max(
                  l_inf_componentwise_residual,
                  residual /
                      (componentwise_residual_offset +
                       std::abs(objective_shard[i] / col_scaling_shard[i])));
            }
          }
        }
        local_dual_correction[shard.Index()] = dual_correction;
        local_dual_full_correction[shard.Index()] = dual_full_correction;
        local_l_inf_residual[shard.Index()] = l_inf_residual;
        local_sumsq_residual[shard.Index()] = sumsq_residual;
        local_l_inf_componentwise_residual[shard.Index()] =
            l_inf_componentwise_residual;
      });
  return ResidualNorms{
      .objective_correction = local_dual_correction.sum(),
      .objective_full_correction = local_dual_full_correction.sum(),
      .l_inf_residual = local_l_inf_residual.lpNorm<Eigen::Infinity>(),
      .l_2_residual = std::sqrt(local_sumsq_residual.sum()),
      .l_inf_componentwise_residual =
          local_l_inf_componentwise_residual.lpNorm<Eigen::Infinity>(),
  };
}

// Returns Qx.
VectorXd ObjectiveProduct(const ShardedQuadraticProgram& sharded_qp,
                          const VectorXd& primal_solution) {
  CHECK_EQ(primal_solution.size(), sharded_qp.PrimalSize());
  VectorXd result(primal_solution.size());
  if (IsLinearProgram(sharded_qp.Qp())) {
    SetZero(sharded_qp.PrimalSharder(), result);
  } else {
    sharded_qp.PrimalSharder().ParallelForEachShard(
        [&](const Sharder::Shard& shard) {
          shard(result) =
              shard(*sharded_qp.Qp().objective_matrix) * shard(primal_solution);
        });
  }
  return result;
}

// Returns 1/2 x^T Q x (the quadratic term in the objective).
double QuadraticObjective(const ShardedQuadraticProgram& sharded_qp,
                          const VectorXd& primal_solution,
                          const VectorXd& objective_product) {
  CHECK_EQ(primal_solution.size(), sharded_qp.PrimalSize());
  CHECK_EQ(objective_product.size(), sharded_qp.PrimalSize());
  return 0.5 *
         Dot(objective_product, primal_solution, sharded_qp.PrimalSharder());
}

// Returns `objective_product` + c − A^T y when `use_zero_primal_objective` is
// false, and returns − A^T y when `use_zero_primal_objective` is true.
// `objective_product` is passed by value, and modified in place.
VectorXd PrimalGradientFromObjectiveProduct(
    const ShardedQuadraticProgram& sharded_qp, const VectorXd& dual_solution,
    VectorXd objective_product, bool use_zero_primal_objective = false) {
  const QuadraticProgram& qp = sharded_qp.Qp();
  CHECK_EQ(dual_solution.size(), sharded_qp.DualSize());
  CHECK_EQ(objective_product.size(), sharded_qp.PrimalSize());

  // Note that this modifies `objective_product`, replacing its entries with
  // the primal gradient.
  sharded_qp.ConstraintMatrixSharder().ParallelForEachShard(
      [&](const Sharder::Shard& shard) {
        if (use_zero_primal_objective) {
          shard(objective_product) =
              -shard(qp.constraint_matrix).transpose() * dual_solution;
        } else {
          shard(objective_product) +=
              shard(qp.objective_vector) -
              shard(qp.constraint_matrix).transpose() * dual_solution;
        }
      });
  return objective_product;
}

// Returns the value of y term in the objective of the dual problem, that is,
// (l^c)^T[y]_+ − (u^c)^T[y]_− in the dual objective from
// https://developers.google.com/optimization/lp/pdlp_math.
double DualObjectiveBoundsTerm(const ShardedQuadraticProgram& sharded_qp,
                               const VectorXd& dual_solution) {
  const QuadraticProgram& qp = sharded_qp.Qp();
  return sharded_qp.DualSharder().ParallelSumOverShards(
      [&](const Sharder::Shard& shard) {
        // This assumes that the dual variables are feasible, that is, that
        // the term corresponding to the "y" variables in the dual objective
        // in https://developers.google.com/optimization/lp/pdlp_math is finite.
        const auto lower_bound_shard = shard(qp.constraint_lower_bounds);
        const auto upper_bound_shard = shard(qp.constraint_upper_bounds);
        const auto dual_shard = shard(dual_solution);
        // Can't use `.dot(.cwiseMin(...))` because that gives 0 * inf = NaN.
        double sum = 0.0;
        for (int64_t i = 0; i < dual_shard.size(); ++i) {
          if (dual_shard[i] > 0.0) {
            sum += lower_bound_shard[i] * dual_shard[i];
          } else if (dual_shard[i] < 0.0) {
            sum += upper_bound_shard[i] * dual_shard[i];
          }
        }
        return sum;
      });
}

// Computes the projection of `vector` onto a pseudo-random vector determined
// by `seed_generator`. `seed_generator` is used as the source of a random seed
// for each shard's portion of the vector.
double RandomProjection(const VectorXd& vector, const Sharder& sharder,
                        std::mt19937& seed_generator) {
  std::vector<std::mt19937> shard_seeds;
  shard_seeds.reserve(sharder.NumShards());
  for (int shard = 0; shard < sharder.NumShards(); ++shard) {
    shard_seeds.emplace_back((seed_generator)());
  }
  // Computes `vector` * gaussian_random_vector and ||gaussian_random_vector||^2
  // to normalize by afterwards.
  VectorXd dot_product(sharder.NumShards());
  VectorXd gaussian_norm_squared(sharder.NumShards());
  sharder.ParallelForEachShard([&](const Sharder::Shard& shard) {
    const auto vector_shard = shard(vector);
    double shard_dot_product = 0.0;
    double shard_norm_squared = 0.0;
    std::mt19937 random{shard_seeds[shard.Index()]};
    for (int64_t i = 0; i < vector_shard.size(); ++i) {
      const double projection_element = absl::Gaussian(random, 0.0, 1.0);
      shard_dot_product += projection_element * vector_shard[i];
      shard_norm_squared += MathUtil::Square(projection_element);
    }
    dot_product[shard.Index()] = shard_dot_product;
    gaussian_norm_squared[shard.Index()] = shard_norm_squared;
  });
  return dot_product.sum() / std::sqrt(gaussian_norm_squared.sum());
}
}  // namespace

ConvergenceInformation ComputeConvergenceInformation(
    const PrimalDualHybridGradientParams& params,
    const ShardedQuadraticProgram& scaled_sharded_qp,
    const Eigen::VectorXd& col_scaling_vec,
    const Eigen::VectorXd& row_scaling_vec,
    const Eigen::VectorXd& scaled_primal_solution,
    const Eigen::VectorXd& scaled_dual_solution,
    const double componentwise_primal_residual_offset,
    const double componentwise_dual_residual_offset, PointType candidate_type) {
  const QuadraticProgram& qp = scaled_sharded_qp.Qp();
  CHECK_EQ(col_scaling_vec.size(), scaled_sharded_qp.PrimalSize());
  CHECK_EQ(row_scaling_vec.size(), scaled_sharded_qp.DualSize());
  CHECK_EQ(scaled_primal_solution.size(), scaled_sharded_qp.PrimalSize());
  CHECK_EQ(scaled_dual_solution.size(), scaled_sharded_qp.DualSize());

  // See https://developers.google.com/optimization/lp/pdlp_math#rescaling for
  // notes describing the connection between the scaled and unscaled problem.

  ConvergenceInformation result;
  ResidualNorms primal_residuals = PrimalResidualNorms(
      scaled_sharded_qp, row_scaling_vec, scaled_primal_solution,
      componentwise_primal_residual_offset);
  result.set_l_inf_primal_residual(primal_residuals.l_inf_residual);
  result.set_l2_primal_residual(primal_residuals.l_2_residual);
  result.set_l_inf_componentwise_primal_residual(
      primal_residuals.l_inf_componentwise_residual);

  result.set_l_inf_primal_variable(
      ScaledLInfNorm(scaled_primal_solution, col_scaling_vec,
                     scaled_sharded_qp.PrimalSharder()));
  result.set_l2_primal_variable(ScaledNorm(scaled_primal_solution,
                                           col_scaling_vec,
                                           scaled_sharded_qp.PrimalSharder()));
  result.set_l_inf_dual_variable(ScaledLInfNorm(
      scaled_dual_solution, row_scaling_vec, scaled_sharded_qp.DualSharder()));
  result.set_l2_dual_variable(ScaledNorm(scaled_dual_solution, row_scaling_vec,
                                         scaled_sharded_qp.DualSharder()));

  VectorXd scaled_objective_product =
      ObjectiveProduct(scaled_sharded_qp, scaled_primal_solution);
  const double quadratic_objective = QuadraticObjective(
      scaled_sharded_qp, scaled_primal_solution, scaled_objective_product);
  VectorXd scaled_primal_gradient = PrimalGradientFromObjectiveProduct(
      scaled_sharded_qp, scaled_dual_solution,
      std::move(scaled_objective_product));
  result.set_primal_objective(qp.ApplyObjectiveScalingAndOffset(
      quadratic_objective + Dot(qp.objective_vector, scaled_primal_solution,
                                scaled_sharded_qp.PrimalSharder())));

  // This is the dual objective from
  // https://developers.google.com/optimization/lp/pdlp_math minus the last term
  // (involving r). All scaling terms cancel out.
  const double dual_objective_piece =
      -quadratic_objective +
      DualObjectiveBoundsTerm(scaled_sharded_qp, scaled_dual_solution);

  ResidualNorms dual_residuals = DualResidualNorms(
      params, scaled_sharded_qp, col_scaling_vec, scaled_primal_solution,
      scaled_primal_gradient, componentwise_dual_residual_offset);
  result.set_dual_objective(qp.ApplyObjectiveScalingAndOffset(
      dual_objective_piece + dual_residuals.objective_correction));
  result.set_corrected_dual_objective(qp.ApplyObjectiveScalingAndOffset(
      dual_objective_piece + dual_residuals.objective_full_correction));
  result.set_l_inf_dual_residual(dual_residuals.l_inf_residual);
  result.set_l2_dual_residual(dual_residuals.l_2_residual);
  result.set_l_inf_componentwise_dual_residual(
      dual_residuals.l_inf_componentwise_residual);

  result.set_candidate_type(candidate_type);
  return result;
}

namespace {

double PrimalRayMaxSignViolation(const ShardedQuadraticProgram& sharded_qp,
                                 const VectorXd& col_scaling_vec,
                                 const VectorXd& scaled_primal_ray) {
  VectorXd primal_ray_local_max_sign_violation(
      sharded_qp.PrimalSharder().NumShards());
  sharded_qp.PrimalSharder().ParallelForEachShard(
      [&](const Sharder::Shard& shard) {
        const auto lower_bound_shard =
            shard(sharded_qp.Qp().variable_lower_bounds);
        const auto upper_bound_shard =
            shard(sharded_qp.Qp().variable_upper_bounds);
        const auto ray_shard = shard(scaled_primal_ray);
        const auto scale_shard = shard(col_scaling_vec);
        double local_max = 0.0;
        for (int64_t i = 0; i < ray_shard.size(); ++i) {
          if (std::isfinite(lower_bound_shard[i])) {
            local_max = std::max(local_max, -ray_shard[i] * scale_shard[i]);
          }
          if (std::isfinite(upper_bound_shard[i])) {
            local_max = std::max(local_max, ray_shard[i] * scale_shard[i]);
          }
        }
        primal_ray_local_max_sign_violation[shard.Index()] = local_max;
      });
  return primal_ray_local_max_sign_violation.lpNorm<Eigen::Infinity>();
}

}  // namespace

InfeasibilityInformation ComputeInfeasibilityInformation(
    const PrimalDualHybridGradientParams& params,
    const ShardedQuadraticProgram& scaled_sharded_qp,
    const Eigen::VectorXd& col_scaling_vec,
    const Eigen::VectorXd& row_scaling_vec,
    const Eigen::VectorXd& scaled_primal_ray,
    const Eigen::VectorXd& scaled_dual_ray,
    const Eigen::VectorXd& primal_solution_for_residual_tests,
    PointType candidate_type) {
  const QuadraticProgram& qp = scaled_sharded_qp.Qp();
  CHECK_EQ(col_scaling_vec.size(), scaled_sharded_qp.PrimalSize());
  CHECK_EQ(row_scaling_vec.size(), scaled_sharded_qp.DualSize());
  CHECK_EQ(scaled_primal_ray.size(), scaled_sharded_qp.PrimalSize());
  CHECK_EQ(scaled_dual_ray.size(), scaled_sharded_qp.DualSize());

  double l_inf_primal = ScaledLInfNorm(scaled_primal_ray, col_scaling_vec,
                                       scaled_sharded_qp.PrimalSharder());
  double l_inf_dual = ScaledLInfNorm(scaled_dual_ray, row_scaling_vec,
                                     scaled_sharded_qp.DualSharder());
  InfeasibilityInformation result;
  // Compute primal infeasibility information.
  VectorXd scaled_primal_gradient = PrimalGradientFromObjectiveProduct(
      scaled_sharded_qp, scaled_dual_ray,
      ZeroVector(scaled_sharded_qp.PrimalSharder()),
      /*use_zero_primal_objective=*/true);
  // We don't use `dual_residuals.l_inf_componentwise_residual`, so don't need
  // to set `componentwise_residual_offset` to a meaningful value.
  ResidualNorms dual_residuals = DualResidualNorms(
      params, scaled_sharded_qp, col_scaling_vec,
      primal_solution_for_residual_tests, scaled_primal_gradient,
      /*componentwise_residual_offset=*/0.0);

  double dual_ray_objective =
      DualObjectiveBoundsTerm(scaled_sharded_qp, scaled_dual_ray) +
      dual_residuals.objective_correction;
  if (l_inf_dual > 0) {
    result.set_dual_ray_objective(dual_ray_objective / l_inf_dual);
    result.set_max_dual_ray_infeasibility(dual_residuals.l_inf_residual /
                                          l_inf_dual);
  } else {
    result.set_dual_ray_objective(0.0);
    result.set_max_dual_ray_infeasibility(0.0);
  }

  // Compute dual infeasibility information. We don't use
  // `primal_residuals.l_inf_componentwise_residual`, so don't need to set
  // `componentwise_residual_offset` to a meaningful value.
  ResidualNorms primal_residuals =
      PrimalResidualNorms(scaled_sharded_qp, row_scaling_vec, scaled_primal_ray,
                          /*componentwise_residual_offset=*/0.0,
                          /*use_homogeneous_constraint_bounds=*/true);

  // The primal ray should have been projected onto the feasibility bounds, so
  // that it has the correct signs.
  DCHECK_EQ(PrimalRayMaxSignViolation(scaled_sharded_qp, col_scaling_vec,
                                      scaled_primal_ray),
            0.0);

  if (l_inf_primal > 0.0) {
    VectorXd scaled_objective_product =
        ObjectiveProduct(scaled_sharded_qp, scaled_primal_ray);
    result.set_primal_ray_quadratic_norm(
        LInfNorm(scaled_objective_product, scaled_sharded_qp.PrimalSharder()) /
        l_inf_primal);
    result.set_max_primal_ray_infeasibility(primal_residuals.l_inf_residual /
                                            l_inf_primal);
    result.set_primal_ray_linear_objective(
        Dot(scaled_primal_ray, qp.objective_vector,
            scaled_sharded_qp.PrimalSharder()) /
        l_inf_primal);
  } else {
    result.set_primal_ray_quadratic_norm(0.0);
    result.set_max_primal_ray_infeasibility(0.0);
    result.set_primal_ray_linear_objective(0.0);
  }

  result.set_candidate_type(candidate_type);
  return result;
}

ConvergenceInformation ComputeScaledConvergenceInformation(
    const PrimalDualHybridGradientParams& params,
    const ShardedQuadraticProgram& sharded_qp, const VectorXd& primal_solution,
    const VectorXd& dual_solution,
    const double componentwise_primal_residual_offset,
    const double componentwise_dual_residual_offset, PointType candidate_type) {
  return ComputeConvergenceInformation(
      params, sharded_qp, OnesVector(sharded_qp.PrimalSharder()),
      OnesVector(sharded_qp.DualSharder()), primal_solution, dual_solution,
      componentwise_primal_residual_offset, componentwise_dual_residual_offset,
      candidate_type);
}

VectorXd ReducedCosts(const PrimalDualHybridGradientParams& params,
                      const ShardedQuadraticProgram& sharded_qp,
                      const VectorXd& primal_solution,
                      const VectorXd& dual_solution,
                      bool use_zero_primal_objective) {
  VectorXd objective_product;
  if (use_zero_primal_objective) {
    objective_product = ZeroVector(sharded_qp.PrimalSharder());
  } else {
    objective_product = ObjectiveProduct(sharded_qp, primal_solution);
  }
  return PrimalGradientFromObjectiveProduct(sharded_qp, dual_solution,
                                            std::move(objective_product),
                                            use_zero_primal_objective);
}

std::optional<ConvergenceInformation> GetConvergenceInformation(
    const IterationStats& stats, PointType candidate_type) {
  for (const auto& convergence_information : stats.convergence_information()) {
    if (convergence_information.candidate_type() == candidate_type) {
      return convergence_information;
    }
  }
  return std::nullopt;
}

std::optional<InfeasibilityInformation> GetInfeasibilityInformation(
    const IterationStats& stats, PointType candidate_type) {
  for (const auto& infeasibility_information :
       stats.infeasibility_information()) {
    if (infeasibility_information.candidate_type() == candidate_type) {
      return infeasibility_information;
    }
  }
  return std::nullopt;
}

std::optional<PointMetadata> GetPointMetadata(const IterationStats& stats,
                                              const PointType point_type) {
  for (const auto& metadata : stats.point_metadata()) {
    if (metadata.point_type() == point_type) {
      return metadata;
    }
  }
  return std::nullopt;
}

void SetRandomProjections(const ShardedQuadraticProgram& sharded_qp,
                          const Eigen::VectorXd& primal_solution,
                          const Eigen::VectorXd& dual_solution,
                          const std::vector<int>& random_projection_seeds,
                          PointMetadata& metadata) {
  for (const int random_projection_seed : random_projection_seeds) {
    std::mt19937 seed_generator(random_projection_seed);
    metadata.mutable_random_primal_projections()->Add(RandomProjection(
        primal_solution, sharded_qp.PrimalSharder(), seed_generator));
    metadata.mutable_random_dual_projections()->Add(RandomProjection(
        dual_solution, sharded_qp.DualSharder(), seed_generator));
  }
}

}  // namespace operations_research::pdlp