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

// optimization_base.hpp
#ifndef SANGI_OPTIMIZATION_BASE_HPP
#define SANGI_OPTIMIZATION_BASE_HPP

#include <math/core/common.hpp>
#include <math/core/traits.hpp>
#include <math/core/vector.hpp>
#include <math/core/matrix.hpp>
#include <math/concepts/algebraic_concepts.hpp>
#include <functional>
#include <limits>
#include <tuple>
#include <cmath>
#include <optional>

namespace sangi {

// 最適化アルゴリズムの共通オプションの構造体
template<concepts::OrderedField T>
struct OptimizationOptions {
    T tolerance = std::numeric_limits<T>::epsilon() * 100;
    size_t max_iterations = 1000;
    bool verbose = false;
};

// 最適化結果の構造体
template<typename PointType, concepts::OrderedField ValueType>
struct OptimizationResult {
    PointType point;       // 最適化された点
    ValueType value;       // 最適化された値
    bool success;          // 最適化の成功フラグ
    size_t iterations;     // 実行された反復回数
    
    // 成功結果を作成する静的ヘルパー
    static OptimizationResult makeSuccess(const PointType& p, const ValueType& v, size_t iters) {
        return {p, v, true, iters};
    }

    // 失敗結果を作成する静的ヘルパー
    static OptimizationResult makeFailure(const PointType& p, const ValueType& v, size_t iters) {
        return {p, v, false, iters};
    }
};

// 一次元問題の結果型
template<concepts::OrderedField T>
using OptimizationResult1D = OptimizationResult<T, T>;

// ライン探索のパラメータ構造体
template<concepts::OrderedField T>
struct LineSearchParams {
    T alpha_init = T(0.1);          // 初期ステップサイズ
    T armijo_const = T(1e-4);       // Armijoルールの定数
    T reduction_factor = T(0.5);    // バックトラックの縮小係数
    size_t max_iterations = 20;     // 最大反復回数
};

// Nelder-Mead 法のオプション
template<concepts::OrderedField T>
struct NelderMeadOptions {
    T tolerance = std::numeric_limits<T>::epsilon() * T(100);
    size_t max_iterations = 1000;
    T initial_step = T(1);          // 初期シンプレックスのステップサイズ
    T reflection = T(1);            // 反射係数 α
    T expansion = T(2);             // 拡大係数 γ
    T contraction = T(0.5);         // 縮小係数 ρ
    T shrink = T(0.5);              // 全体縮小係数 σ
};

// Gauss-Newton 法のオプション
template<concepts::OrderedField T>
struct GaussNewtonOptions {
    T tolerance = std::numeric_limits<T>::epsilon() * T(100);
    size_t max_iterations = 100;
    T damping = T(1e-6);            // LM 正則化パラメータ λ
    bool use_line_search = true;    // Armijo ライン探索を使用するか
};

// 共役勾配法のオプション
template<concepts::OrderedField T>
struct ConjugateGradientOptions {
    T tolerance = std::numeric_limits<T>::epsilon() * T(100);
    size_t max_iterations = 1000;
    size_t restart_interval = 0;    // 0 = n（次元数）でリスタート
    bool use_polak_ribiere = true;  // true: Polak-Ribière, false: Fletcher-Reeves
};

// L-BFGS のオプション
template<concepts::OrderedField T>
struct LBFGSOptions {
    T tolerance = std::numeric_limits<T>::epsilon() * T(100);
    size_t max_iterations = 1000;
    size_t memory_size = 10;         // 保持する過去ベクトル数 m
};

// L-BFGS-B のオプション（境界制約付き）
template<concepts::OrderedField T>
struct LBFGSBOptions {
    T tolerance = std::numeric_limits<T>::epsilon() * T(100);
    size_t max_iterations = 1000;
    size_t memory_size = 10;
};

// Hooke-Jeeves パターン探索のオプション
template<concepts::OrderedField T>
struct HookeJeevesOptions {
    T tolerance = T(1e-12);          // 評価値変化の収束判定
    size_t max_iterations = 10000;   // 最大評価回数
    T initial_step = T(1);           // 初期ステップサイズ
    T min_step = T(1e-10);           // ステップがこれ以下で収束
};

// 制約付き最適化の共通オプション
template<concepts::OrderedField T>
struct ConstrainedOptions {
    T tolerance = std::numeric_limits<T>::epsilon() * T(100);
    size_t max_iterations = 100;         // 外側ループ
    size_t inner_max_iterations = 1000;  // 内側ソルバー
    T initial_penalty = T(1);            // 初期ペナルティパラメータ μ
    T penalty_growth = T(10);            // μ の成長率
};

// 遺伝的アルゴリズム (実数コード型) のオプション
template<concepts::OrderedField T>
struct GeneticAlgorithmOptions {
    T tolerance = T(1e-8);           // 最良個体の変化が tolerance 以下で収束
    size_t max_generations = 1000;   // 最大世代数
    size_t population_size = 50;     // 個体数
    T crossover_rate = T(0.8);       // 交叉率
    T mutation_rate = T(0.1);        // 突然変異率
    T mutation_scale = T(0.1);       // 突然変異の標準偏差（探索範囲に対する比率）
    size_t tournament_size = 3;      // トーナメント選択のサイズ
    size_t elite_count = 2;          // エリート保存数
    unsigned int seed = 42;          // 乱数シード
};

// 焼きなまし法 (Simulated Annealing) のオプション
template<concepts::OrderedField T>
struct SimulatedAnnealingOptions {
    T tolerance = T(1e-10);              // 最良値変化の収束判定
    size_t max_iterations = 100000;      // 最大反復回数
    T initial_temperature = T(100);      // 初期温度
    T cooling_rate = T(0.995);           // 冷却率: T(k+1) = cooling_rate * T(k)
    T neighbor_scale = T(0.1);           // 近傍生成のガウスノイズ標準偏差（範囲比）
    unsigned int seed = 42;              // 乱数シード
};

// 差分進化法 (Differential Evolution) のオプション
template<concepts::OrderedField T>
struct DifferentialEvolutionOptions {
    T tolerance = T(1e-8);               // 最良個体の変化が tolerance 以下で収束
    size_t max_generations = 1000;       // 最大世代数
    size_t population_size = 50;         // 個体数
    T mutation_factor = T(0.8);          // 差分重み F
    T crossover_rate = T(0.9);           // 交叉確率 CR
    unsigned int seed = 42;              // 乱数シード
};

// 非線形最小二乗フィッティング (Levenberg-Marquardt) のオプション
template<concepts::OrderedField T>
struct LeastSquaresOptions {
    T ftol = T(1e-10);              // 残差二乗和の相対変化
    T xtol = T(1e-10);             // パラメータの相対変化
    T gtol = T(1e-10);             // 勾配ノルムの収束判定
    std::size_t maxIterations = 200;
    T lambdaInit = T(1e-3);        // 初期 LM パラメータ λ
    T lambdaUp = T(10);            // λ 増大率 (ステップ棄却時)
    T lambdaDown = T(0.1);         // λ 縮小率 (ステップ採択時)
};

// 非線形最小二乗フィッティングの結果 (統計量付き)
template<concepts::OrderedField T>
struct LeastSquaresFitResult {
    Vector<T> parameters;          // 最適パラメータ
    Vector<T> residuals;           // 最終残差ベクトル
    Matrix<T> covariance;          // パラメータ共分散行列 s² (J^T J)^{-1}
    Vector<T> standardErrors;      // 各パラメータの標準誤差
    T chiSquared;                  // 残差二乗和 Σ r_i²
    T reducedChiSquared;           // χ² / (m - n)
    std::size_t iterations;        // 反復回数
    bool success;                  // 収束フラグ
    std::size_t dataPoints;        // データ点数 m
    std::size_t numParameters;     // パラメータ数 n
};

// CMA-ES (共分散行列適応進化戦略) のオプション
template<concepts::OrderedField T>
struct CMAESOptions {
    T tolerance = T(1e-10);              // 関数値変化の収束判定 (TolFun)
    T tolx = T(1e-12);                   // 平均ベクトル変化の収束判定 (TolX)
    size_t max_generations = 10000;      // 最大世代数
    size_t population_size = 0;          // λ (0 = 自動: 4 + floor(3*ln(n)))
    T initial_sigma = T(0.3);            // 初期ステップサイズ σ (探索範囲に対する比率)
    T condition_limit = T(1e14);         // 共分散行列の条件数上限 (ConditionCov)
    unsigned int seed = 42;              // 乱数シード
};

// jSO (改良版適応的差分進化法) のオプション
template<concepts::OrderedField T>
struct JSOOptions {
    T tolerance = T(1e-8);               // 最良個体の変化が tolerance 以下で収束
    size_t max_evaluations = 0;          // 最大関数評価数 (0 = 自動: 10000*D)
    size_t population_size = 0;          // NP_init (0 = 自動: 25*D)
    size_t min_population_size = 4;      // NP_min (最小個体数)
    size_t history_size = 5;             // 成功パラメータ記憶の長さ H
    T archive_rate = T(2.6);             // アーカイブサイズ倍率 (NP_init * archive_rate)
    T p_best_rate = T(0.11);             // pbest 選択の上位割合
    unsigned int seed = 42;              // 乱数シード
};

// MKL利用可能性チェック
inline constexpr bool has_mkl_support = SANGI_HAS_MKL == 1;

} // namespace sangi

#endif // SANGI_OPTIMIZATION_BASE_HPP