// Copyright (C) 2026 Kiyotsugu Arai
// SPDX-License-Identifier: LGPL-3.0-or-later

// linear_programming_impl.hpp — implementation of simplex, simplexMaximize
//
// ★2026-06: the old two-phase simplex (detail::SimplexSolver) was confirmed by a codex audit to have a *correctness bug*
//   (Phase 1 erroneously dropped rows still containing artificial variables, losing real constraints → returns a constraint-violating solution with status=Optimal).
//   Keeping the public API (LinearProgram→LPResult) intact, the internals were replaced with delegation to the robust robust_lp (Mehrotra interior-point method,
//   with finite guards, "the returned solution always satisfies Ax=b, x≥0").
//   ★Behavioral differences (callers beware): being an interior-point method, robust_lp (a) returns a non-basic *interior* optimal solution
//   (not a vertex/basic solution; the analytic center on non-unique optimal faces), (b) has accuracy roughly tol~1e-7 (not the exact-vertex 1e-12 level of the old simplex),
//   (c) computes in double (T is internally converted to double). See the top of robust_lp.hpp for details.

#ifndef SANGI_LINEAR_PROGRAMMING_IMPL_HPP
#define SANGI_LINEAR_PROGRAMMING_IMPL_HPP

#include <math/optimization/linear_programming.hpp>
#include <math/optimization/robust_lp.hpp>

#include <vector>

namespace sangi {

template<typename T>
LPResult<T> simplex(
    const LinearProgram<T>& lp,
    T eps,
    size_t maxIterations)
{
    (void)eps;
    const size_t nvar = lp.numVariables();
    const size_t mcon = lp.numConstraints();
    if (nvar == 0 || mcon == 0) return LPResult<T>::infeasible(0);

    // Convert LinearProgram (≤/≥/=, x≥0) to the robust_lp standard form min c·x s.t. Ax=b, x≥0.
    //   ≤: +slack, ≥: −slack, =: no slack. x=[original variables nvar | slack nslack].
    size_t nslack = 0;
    for (const auto& con : lp.constraints)
        if (con.type != ConstraintType::Equal) ++nslack;
    const size_t N = nvar + nslack;

    std::vector<std::vector<double>> A(mcon, std::vector<double>(N, 0.0));
    std::vector<double> b(mcon, 0.0), c(N, 0.0);
    for (size_t j = 0; j < nvar; ++j) c[j] = static_cast<double>(lp.objective[j]);
    size_t sidx = 0;
    for (size_t i = 0; i < mcon; ++i) {
        const auto& con = lp.constraints[i];
        for (size_t j = 0; j < nvar && j < con.coefficients.size(); ++j)
            A[i][j] = static_cast<double>(con.coefficients[j]);
        if (con.type == ConstraintType::LessEqual)      A[i][nvar + sidx++] =  1.0;
        else if (con.type == ConstraintType::GreaterEqual) A[i][nvar + sidx++] = -1.0;
        b[i] = static_cast<double>(con.rhs);
    }

    rlp::Options opt;
    opt.tol = 1e-10;                                   // tight, aiming for vertex-level accuracy
    opt.maxIter = (maxIterations > 0 && maxIterations < 400) ? static_cast<int>(maxIterations) : 200;
    auto res = rlp::solveStandardLP(c, A, b, opt);

    const T objConst = lp.objectiveConstant;
    const T objVal = static_cast<T>(res.objective) + objConst;
    std::vector<T> sol(nvar);
    for (size_t j = 0; j < nvar && j < res.x.size(); ++j) sol[j] = static_cast<T>(res.x[j]);

    switch (res.status) {
        case rlp::Status::Optimal:   return LPResult<T>::optimal(objVal, sol, res.iterations);
        case rlp::Status::Unbounded: return LPResult<T>::unbounded(res.iterations);
        case rlp::Status::Infeasible:return LPResult<T>::infeasible(res.iterations);
        default:                     return LPResult<T>::maxIterations(objVal, sol, res.iterations);
    }
}

template<typename T>
LPResult<T> simplexMaximize(
    const LinearProgram<T>& lp,
    T eps,
    size_t maxIterations)
{
    // copy with the objective function sign negated
    LinearProgram<T> minLP = lp;
    for (auto& c : minLP.objective) c = -c;
    minLP.objectiveConstant = -lp.objectiveConstant;

    auto result = simplex(minLP, eps, maxIterations);
    result.objectiveValue = -result.objectiveValue;   // restore the sign of the objective value
    return result;
}

} // namespace sangi

#endif // SANGI_LINEAR_PROGRAMMING_IMPL_HPP