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

// linear_programming.hpp
// 線形計画法 (シンプレックス法)
//
// 2フェーズ改訂シンプレックス法による線形計画問題の求解
// 参考: 奥村「C言語によるアルゴリズム事典」simplex.c を元に
//       型安全なテンプレート実装に全面改良
//
// 改良点:
//   - グローバル変数 → 構造体にカプセル化
//   - 固定サイズ配列 → 動的 Matrix<T>/Vector<T>
//   - scanf → 関数引数として問題を受け取る
//   - エラー処理: exit() → 結果構造体で状態を返す
//   - 等式制約 (=), 不等式制約 (<=, >=) すべて対応

#ifndef SANGI_LINEAR_PROGRAMMING_HPP
#define SANGI_LINEAR_PROGRAMMING_HPP

#include "optimization_base.hpp"
#include <math/core/vector.hpp>
#include <math/core/matrix.hpp>

#include <vector>
#include <cmath>
#include <algorithm>
#include <limits>
#include <cassert>
#include <stdexcept>

namespace sangi {

/// 制約の不等号の種類
enum class ConstraintType {
    LessEqual,    // <=
    GreaterEqual, // >=
    Equal         // =
};

/// 線形計画問題の制約条件 1行分
template<typename T>
struct LinearConstraint {
    std::vector<T> coefficients;  // 左辺の係数 (変数の数と同じ長さ)
    ConstraintType type;          // 不等号の種類
    T rhs;                        // 右辺定数

    LinearConstraint() : type(ConstraintType::LessEqual), rhs(T(0)) {}

    LinearConstraint(std::vector<T> coeffs, ConstraintType t, T r)
        : coefficients(std::move(coeffs)), type(t), rhs(r) {}
};

/// 線形計画問題の定義
/// minimize c^T x   subject to  A_i x {<=, >=, =} b_i
template<typename T>
struct LinearProgram {
    std::vector<T> objective;                      // 目的関数の係数 (最小化)
    T objectiveConstant = T(0);                    // 目的関数の定数項
    std::vector<LinearConstraint<T>> constraints;  // 制約条件のリスト

    size_t numVariables() const { return objective.size(); }
    size_t numConstraints() const { return constraints.size(); }
};

/// LP ソルバーの状態
enum class LPStatus {
    Optimal,       // 最適解が見つかった
    Infeasible,    // 実行可能解なし
    Unbounded,     // 目的関数が非有界
    MaxIterations  // 最大反復回数に到達
};

/// LP ソルバーの結果
template<typename T>
struct LPResult {
    LPStatus status;
    T objectiveValue;              // 最適な目的関数値
    std::vector<T> solution;       // 最適解 (各変数の値)
    size_t iterations;             // 総ピボット回数

    static LPResult optimal(T val, std::vector<T> sol, size_t iters) {
        return {LPStatus::Optimal, val, std::move(sol), iters};
    }
    static LPResult infeasible(size_t iters) {
        return {LPStatus::Infeasible, T(0), {}, iters};
    }
    static LPResult unbounded(size_t iters) {
        return {LPStatus::Unbounded, -std::numeric_limits<T>::infinity(), {}, iters};
    }
    static LPResult maxIterations(T val, std::vector<T> sol, size_t iters) {
        return {LPStatus::MaxIterations, val, std::move(sol), iters};
    }
};

namespace detail {

/// 2フェーズ改訂シンプレックス法の内部実装
template<typename T>
class SimplexSolver {
public:
    SimplexSolver(const LinearProgram<T>& lp, T eps, size_t maxIter)
        : eps_(eps), maxIter_(maxIter), totalPivots_(0)
    {
        setup(lp);
    }

    LPResult<T> solve() {
        // フェーズ 1: 人工変数がある場合、実行可能基底解を見つける
        if (n3_ > n2_) {
            auto status = phase1();
            if (status == LPStatus::Infeasible) {
                return LPResult<T>::infeasible(totalPivots_);
            }
        }

        // フェーズ 2: 目的関数を最適化する
        auto status = phase2();
        if (status == LPStatus::Unbounded) {
            return LPResult<T>::unbounded(totalPivots_);
        }
        if (status == LPStatus::MaxIterations) {
            return LPResult<T>::maxIterations(
                extractObjectiveValue(), extractSolution(), totalPivots_);
        }

        return LPResult<T>::optimal(
            extractObjectiveValue(), extractSolution(), totalPivots_);
    }

private:
    T eps_;
    size_t maxIter_;
    size_t totalPivots_;

    size_t m_;    // 制約の数
    size_t n_;    // 元の変数の数
    size_t n1_;   // n + >= スラック変数の数
    size_t n2_;   // n1 + <= スラック変数の数
    size_t n3_;   // n2 + 人工変数の数
    size_t jmax_; // タブローの最右列番号

    // 制約係数行列 a[0..m][0..n] — 行0は目的関数
    std::vector<std::vector<T>> a_;

    // 変換行列 q[0..m][0..m]
    std::vector<std::vector<T>> q_;

    // ピボット列の一時バッファ
    std::vector<T> pivotCol_;

    // col[i] = 行 i の基本変数の列番号 (i=1..m)
    std::vector<int> col_;

    // row[j] = 列 j が基本変数なら対応する行番号、非基本なら 0
    std::vector<int> row_;

    // nonzeroRow[j] = スラック/人工変数 j が非零の行番号
    std::vector<int> nonzeroRow_;

    // 各制約の不等号種別
    std::vector<ConstraintType> ineq_;

    // 目的関数の係数 c[0..n] (c[0]=-定数項, c[1..n]=各変数の係数)
    // フェーズ2 開始時に a_[0] にコピーする (元のCコードと同じ構造)
    std::vector<T> c_;

    // 目的関数の定数項
    T objConst_;

    /// 問題を内部形式にセットアップ
    void setup(const LinearProgram<T>& lp) {
        m_ = lp.numConstraints();
        n_ = lp.numVariables();
        objConst_ = lp.objectiveConstant;

        // 右辺を非負にする前処理
        ineq_.resize(m_ + 1);
        a_.assign(m_ + 1, std::vector<T>(n_ + 1, T(0)));

        // 目的関数の係数を保存 (a_[0] はフェーズ2で設定)
        c_.assign(n_ + 1, T(0));
        for (size_t j = 0; j < n_; ++j) {
            c_[j + 1] = lp.objective[j];
        }
        c_[0] = -objConst_;

        // 制約条件を格納 (右辺が負の場合は不等号を反転)
        for (size_t i = 0; i < m_; ++i) {
            const auto& c = lp.constraints[i];
            assert(c.coefficients.size() == n_);

            ineq_[i + 1] = c.type;
            a_[i + 1][0] = c.rhs;

            for (size_t j = 0; j < n_; ++j) {
                a_[i + 1][j + 1] = c.coefficients[j];
            }

            // 右辺を非負に正規化
            if (a_[i + 1][0] < T(0)) {
                flipInequality(i + 1);
                for (size_t j = 0; j <= n_; ++j) {
                    a_[i + 1][j] = -a_[i + 1][j];
                }
            } else if (a_[i + 1][0] == T(0) && ineq_[i + 1] == ConstraintType::GreaterEqual) {
                // 0 >= ... は <= に反転 (係数のみ反転、右辺は0のまま)
                ineq_[i + 1] = ConstraintType::LessEqual;
                for (size_t j = 1; j <= n_; ++j) {
                    a_[i + 1][j] = -a_[i + 1][j];
                }
            }
        }

        // スラック変数と人工変数の準備
        prepare();
    }

    /// 不等号を反転する
    void flipInequality(size_t i) {
        if (ineq_[i] == ConstraintType::LessEqual) {
            ineq_[i] = ConstraintType::GreaterEqual;
        } else if (ineq_[i] == ConstraintType::GreaterEqual) {
            ineq_[i] = ConstraintType::LessEqual;
        }
        // Equal はそのまま
    }

    /// スラック変数と人工変数を導入する
    void prepare() {
        size_t totalCols = n_ + 2 * m_ + 1;

        col_.assign(m_ + 1, 0);
        row_.assign(totalCols, 0);
        nonzeroRow_.assign(totalCols, 0);

        // >= 制約に -1 のスラック変数を追加
        n1_ = n_;
        for (size_t i = 1; i <= m_; ++i) {
            if (ineq_[i] == ConstraintType::GreaterEqual) {
                ++n1_;
                nonzeroRow_[n1_] = static_cast<int>(i);
            }
        }

        // <= 制約に +1 のスラック変数を追加 (初期基本変数として使用)
        n2_ = n1_;
        for (size_t i = 1; i <= m_; ++i) {
            if (ineq_[i] == ConstraintType::LessEqual) {
                ++n2_;
                col_[i] = static_cast<int>(n2_);
                nonzeroRow_[n2_] = static_cast<int>(i);
                row_[n2_] = static_cast<int>(i);
            }
        }

        // >= と = の制約に人工変数を追加
        n3_ = n2_;
        for (size_t i = 1; i <= m_; ++i) {
            if (ineq_[i] != ConstraintType::LessEqual) {
                ++n3_;
                col_[i] = static_cast<int>(n3_);
                nonzeroRow_[n3_] = static_cast<int>(i);
                row_[n3_] = static_cast<int>(i);
            }
        }

        // 変換行列 Q を単位行列に初期化
        q_.assign(m_ + 1, std::vector<T>(m_ + 1, T(0)));
        for (size_t i = 0; i <= m_; ++i) {
            q_[i][i] = T(1);
        }

        pivotCol_.resize(m_ + 1);
    }

    /// タブロー要素 (i, j) を計算する
    /// 改訂シンプレックス法: 完全なタブローを保持せず、Q*A で動的に計算
    T tableau(size_t i, size_t j) const {
        if (col_[i] < 0) return T(0); // 削除された行

        if (j <= n_) {
            // 元の変数列: Q の i 行と A の j 列の内積
            T s = T(0);
            for (size_t k = 0; k <= m_; ++k) {
                s += q_[i][k] * a_[k][j];
            }
            return s;
        }

        // スラック/人工変数列
        T s = q_[i][nonzeroRow_[j]];
        if (j <= n1_) return -s;       // >= のスラック変数 (係数 -1)
        if (j <= n2_ || i != 0) return s;  // <= のスラック変数 (係数 +1)
        return s + T(1);               // 人工変数 (目的関数行に +1 補正)
    }

    /// ピボット操作: 行 ip の基本変数を列 jp の変数に置き換える
    /// 注: Q 行列の列 0 は操作しない (目的関数の定数項を保護)
    void pivot(size_t ip, size_t jp) {
        T u = pivotCol_[ip];

        // ピボット行を正規化 (列 1..m のみ)
        for (size_t j = 1; j <= m_; ++j) {
            q_[ip][j] /= u;
        }

        // 他の行からピボット行の寄与を除去 (列 1..m のみ)
        for (size_t i = 0; i <= m_; ++i) {
            if (i != ip) {
                u = pivotCol_[i];
                for (size_t j = 1; j <= m_; ++j) {
                    q_[i][j] -= q_[ip][j] * u;
                }
            }
        }

        // 基本変数の索引を更新
        row_[col_[ip]] = 0;
        col_[ip] = static_cast<int>(jp);
        row_[jp] = static_cast<int>(ip);
    }

    /// 最小化ループ: 非基本変数から被約費用が負のものを選んでピボット
    LPStatus minimize() {
        for (;;) {
            if (totalPivots_ >= maxIter_) {
                return LPStatus::MaxIterations;
            }

            // ピボット列の選択: 被約費用が最も負の非基本変数
            size_t jp = 0;
            T bestCost = -eps_;
            for (size_t j = 1; j <= jmax_; ++j) {
                if (row_[j] == 0) {
                    T cost = tableau(0, j);
                    if (cost < bestCost) {
                        bestCost = cost;
                        jp = j;
                    }
                }
            }
            if (jp == 0) break; // すべての被約費用が非負 → 最適

            // ピボット列を計算してピボット行を選択 (最小比テスト)
            for (size_t i = 0; i <= m_; ++i) {
                pivotCol_[i] = tableau(i, jp);
            }

            size_t ip = 0;
            T minRatio = std::numeric_limits<T>::max();
            for (size_t i = 1; i <= m_; ++i) {
                if (pivotCol_[i] > eps_) {
                    T ratio = tableau(i, 0) / pivotCol_[i];
                    if (ratio < minRatio) {
                        minRatio = ratio;
                        ip = i;
                    }
                }
            }

            if (ip == 0) {
                // ピボット列の全要素が非正 → 目的関数は非有界
                return LPStatus::Unbounded;
            }

            pivot(ip, jp);
            ++totalPivots_;
        }

        return LPStatus::Optimal;
    }

    /// フェーズ 1: 人工変数を基底から追い出し、実行可能基底解を見つける
    LPStatus phase1() {
        jmax_ = n3_;

        // 人工変数の和を最小化する目的関数を設定
        // 行 0 の Q 行を修正: 人工変数を持つ行に -1 を加算
        for (size_t i = 0; i <= m_; ++i) {
            if (col_[i] > static_cast<int>(n2_)) {
                q_[0][i] -= T(1);
            }
        }

        auto status = minimize();
        if (status == LPStatus::Unbounded) {
            return LPStatus::Infeasible;
        }

        // 人工変数の目的関数値が正なら実行不能
        T artObj = tableau(0, 0);
        if (artObj < -eps_) {
            return LPStatus::Infeasible;
        }

        // 基底に残っている人工変数の行を冗長として削除
        for (size_t i = 1; i <= m_; ++i) {
            if (col_[i] > static_cast<int>(n2_)) {
                col_[i] = -1; // 行を削除マーク
            }
        }

        // 元の目的関数に復元
        q_[0].assign(m_ + 1, T(0));
        q_[0][0] = T(1);

        // 基底に入っている元変数の被約費用を相殺
        for (size_t i = 1; i <= m_; ++i) {
            int j = col_[i];
            if (j > 0 && j <= static_cast<int>(n_)) {
                T u = c_[j]; // 目的関数の j 列の係数
                if (u != T(0)) {
                    for (size_t k = 1; k <= m_; ++k) {
                        q_[0][k] -= q_[i][k] * u;
                    }
                }
            }
        }

        return LPStatus::Optimal;
    }

    /// フェーズ 2: 元の目的関数を最小化する
    LPStatus phase2() {
        jmax_ = n2_; // 人工変数列は除外
        // 目的関数の係数を a_[0] にコピー (元のCコードの phase2() と同じ)
        for (size_t j = 0; j <= n_; ++j) {
            a_[0][j] = c_[j];
        }
        return minimize();
    }

    /// 最適な目的関数値を取得
    T extractObjectiveValue() const {
        return -tableau(0, 0) + objConst_;
    }

    /// 最適解の変数値を取得
    std::vector<T> extractSolution() const {
        std::vector<T> sol(n_, T(0));
        for (size_t j = 1; j <= n_; ++j) {
            int i = row_[j];
            if (i != 0) {
                sol[j - 1] = tableau(i, 0);
            }
        }
        return sol;
    }
};

} // namespace detail

/// シンプレックス法で線形計画問題を解く
///
/// @tparam T 数値型 (double, float など)
/// @param lp  線形計画問題の定義 (最小化問題)
/// @param eps 許容誤差 (ピボット選択の閾値)
/// @param maxIterations 最大ピボット回数
/// @return 結果 (状態、目的関数値、解ベクトル、反復回数)
///
/// 使用例:
/// @code
///   // minimize -x - 2y  subject to  x + y <= 4,  x <= 3,  y <= 3
///   LinearProgram<double> lp;
///   lp.objective = {-1.0, -2.0};
///   lp.constraints = {
///       {{1.0, 1.0}, ConstraintType::LessEqual, 4.0},
///       {{1.0, 0.0}, ConstraintType::LessEqual, 3.0},
///       {{0.0, 1.0}, ConstraintType::LessEqual, 3.0}
///   };
///   auto result = simplex(lp);
///   // result.status == LPStatus::Optimal
///   // result.objectiveValue == -7.0
///   // result.solution == {1.0, 3.0}
/// @endcode
template<typename T>
LPResult<T> simplex(
    const LinearProgram<T>& lp,
    T eps = T(1e-10),
    size_t maxIterations = 10000);

/// 最大化問題用ラッパー: maximize c^T x → minimize -c^T x
template<typename T>
LPResult<T> simplexMaximize(
    const LinearProgram<T>& lp,
    T eps = T(1e-10),
    size_t maxIterations = 10000);

// 実装は linear_programming_impl.hpp に分離。

} // namespace sangi

#endif // SANGI_LINEAR_PROGRAMMING_HPP