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

// optimization_nd.hpp
#ifndef SANGI_OPTIMIZATION_ND_HPP
#define SANGI_OPTIMIZATION_ND_HPP

#include "optimization_base.hpp"
#include <math/core/vector.hpp>
#include <math/core/matrix.hpp>
#include <math/linalg/solvers.hpp>
#include <math/linalg/decomposition.hpp>
#include <random>
#include <algorithm>
#include <numbers>

namespace sangi {

/**
 * @brief バックトラッキングライン探索
 * @tparam T 順序構造を持つ体型（実数など）
 * @tparam V ベクトル型
 * @param f 最小化する関数
 * @param x 現在の点
 * @param direction 探索方向
 * @param fx 現在の点での関数値
 * @param gradient 現在の点での勾配
 * @param params ライン探索パラメータ
 * @return (アルファ, 新しい点, 新しい関数値, 成功フラグ)のタプル
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
std::tuple<T, V, T, bool> backtracking_line_search(
    const std::function<T(const V&)>& f,
    const V& x,
    const V& direction,
    T fx,
    const V& gradient,
    const LineSearchParams<T>& params = LineSearchParams<T>());

/**
 * @brief ラインサーチ（最急降下法）
 * @tparam T 順序構造を持つ体型（実数など）
 * @tparam V ベクトル型
 * @param f 最小化する関数
 * @param df 関数の勾配を計算する関数
 * @param x0 初期点
 * @param options 最適化オプション
 * @param line_search_params ライン探索パラメータ
 * @return 最小値を与える点と最小値、および成功フラグのタプル
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> steepest_descent(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const V& x0,
    const OptimizationOptions<T>& options = OptimizationOptions<T>(),
    const LineSearchParams<T>& line_search_params = LineSearchParams<T>());

/**
 * @brief 準ニュートン法（BFGS法）
 * @tparam T 順序構造を持つ体型（実数など）
 * @tparam V ベクトル型
 * @tparam M 行列型
 * @param f 最小化する関数
 * @param df 関数の勾配を計算する関数
 * @param x0 初期点
 * @param options 最適化オプション
 * @param line_search_params ライン探索パラメータ
 * @return 最小値を与える点と最小値、および成功フラグのタプル
 */
template<concepts::OrderedField T, typename V, typename M>
    requires concepts::VectorOf<V, T> && concepts::MatrixOf<M, T>
OptimizationResult<V, T> bfgs_minimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const V& x0,
    const OptimizationOptions<T>& options = OptimizationOptions<T>(),
    const LineSearchParams<T>& line_search_params = LineSearchParams<T>());

/**
 * @brief ニュートン共役勾配法（Truncated Newton Method）
 * @tparam T 順序構造を持つ体型（実数など）
 * @tparam V ベクトル型
 * @param f 最小化する関数
 * @param df 関数の勾配を計算する関数
 * @param d2f 関数のヘッセ行列と方向ベクトルの積を計算する関数
 * @param x0 初期点
 * @param options 最適化オプション
 * @param line_search_params ライン探索パラメータ
 * @return 最小値を与える点と最小値、および成功フラグのタプル
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> newton_cg_minimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const std::function<V(const V&, const V&)>& d2f,
    const V& x0,
    const OptimizationOptions<T>& options = OptimizationOptions<T>(),
    const LineSearchParams<T>& line_search_params = LineSearchParams<T>());

/**
 * @brief ベクトルの外積（結果は行列）
 * @tparam T 体型（実数など）
 * @tparam V ベクトル型
 * @tparam M 行列型
 * @param a ベクトル
 * @param b ベクトル
 * @return 外積 a⊗b (要素 a_i * b_j の行列)
 */
template<concepts::Field T, typename V, typename M>
    requires concepts::VectorOf<V, T> && concepts::MatrixOf<M, T>
M outer_product(const V& a, const V& b);

// =====================================================================
// Nelder-Mead シンプレックス法（微分不要の最適化）
// =====================================================================

/**
 * @brief Nelder-Mead シンプレックス法
 * @tparam T 順序体型
 * @tparam V ベクトル型
 * @param f 最小化する関数
 * @param x0 初期点
 * @param options Nelder-Mead オプション
 * @return 最適化結果
 *
 * @note lib++20 からの改善:
 *   - 収束判定: 関数値スプレッド AND シンプレックス直径の二重判定
 *   - タイムアウト時: 最良点を返却（lib++20 は TYPE(0) を返却）
 *   - 標準的な Nelder-Mead 係数 (α=1, γ=2, ρ=0.5, σ=0.5)
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> nelderMeadMinimize(
    const std::function<T(const V&)>& f,
    const V& x0,
    const NelderMeadOptions<T>& options = NelderMeadOptions<T>());

// =====================================================================
// Gauss-Newton 法（非線形最小二乗）
// =====================================================================

/**
 * @brief Gauss-Newton 法（Levenberg-Marquardt 正則化付き）
 * @tparam T 順序体型
 * @tparam V ベクトル型
 * @tparam M 行列型
 * @param residuals 残差関数 r(x): V→V
 * @param jacobian ヤコビ行列関数 J(x): V→M
 * @param x0 初期パラメータ
 * @param options Gauss-Newton オプション
 * @param line_search_params ライン探索パラメータ
 * @return 最適化結果（value は残差の二乗和）
 *
 * @note lib++20 からの改善:
 *   - LM 減衰 (J^T J + λI) で特異性回避
 *   - 古い dX へのフォールバック廃止 → 失敗として報告
 *   - Armijo ライン探索追加（固定ステップではなくなった）
 *   - stdout デバッグ出力を削除
 */
template<concepts::OrderedField T, typename V, typename M>
    requires concepts::VectorOf<V, T> && concepts::MatrixOf<M, T>
OptimizationResult<V, T> gaussNewtonMinimize(
    const std::function<V(const V&)>& residuals,
    const std::function<M(const V&)>& jacobian,
    const V& x0,
    const GaussNewtonOptions<T>& options = GaussNewtonOptions<T>(),
    const LineSearchParams<T>& line_search_params = LineSearchParams<T>());

// =====================================================================
// 非線形共役勾配法（Fletcher-Reeves / Polak-Ribière）
// =====================================================================

/**
 * @brief 非線形共役勾配法
 * @tparam T 順序体型
 * @tparam V ベクトル型
 * @param f 最小化する関数
 * @param df 勾配関数
 * @param x0 初期点
 * @param options 共役勾配法オプション
 * @param line_search_params ライン探索パラメータ
 * @return 最適化結果
 *
 * @note lib++20 からの改善:
 *   - Armijo backtracking（lib++20 は α=0.000001 固定 + 無制限ブラケット拡張）
 *   - n ステップごとのリスタート（lib++20 は共役性喪失を放置）
 *   - β 計算のゼロ除算ガード
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
OptimizationResult<V, T> conjugateGradientMinimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const V& x0,
    const ConjugateGradientOptions<T>& options = ConjugateGradientOptions<T>(),
    const LineSearchParams<T>& line_search_params = LineSearchParams<T>());

// =====================================================================
// L-BFGS (Limited-memory BFGS)
// =====================================================================

/**
 * @brief L-BFGS 法（省メモリ準ニュートン法）
 *
 * BFGS の省メモリ版。ヘッセ逆行列を明示的に保持せず、
 * 過去 m 個の勾配差分で二重ループ再帰により H⁻¹g を O(mn) で計算。
 * 大規模問題（変数 1000+）で標準的に使われる。
 *
 * @param f 最小化する関数
 * @param df 勾配関数
 * @param x0 初期点
 * @param options L-BFGS オプション (memory_size = 保持するベクトル数)
 * @param line_search_params ライン探索パラメータ
 * @return 最適化結果
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> lbfgsMinimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const V& x0,
    const LBFGSOptions<T>& options = LBFGSOptions<T>(),
    const LineSearchParams<T>& line_search_params = LineSearchParams<T>());

// =====================================================================
// L-BFGS-B（境界制約付き L-BFGS）
// =====================================================================

/**
 * @brief L-BFGS-B 法（境界制約付き省メモリ準ニュートン法）
 *
 * 各変数に上下限制約 l_i ≤ x_i ≤ u_i を持つ最適化。
 * 勾配射影法で探索方向を実行可能領域に投影し、
 * 自由変数のみで L-BFGS ステップを計算。
 *
 * @param f 最小化する関数
 * @param df 勾配関数
 * @param x0 初期点（実行可能でなければ自動投影）
 * @param lower 下限（-inf の場合は -numeric_limits::max()）
 * @param upper 上限（+inf の場合は +numeric_limits::max()）
 * @param options L-BFGS-B オプション
 * @param line_search_params ライン探索パラメータ
 * @return 最適化結果
 *
 * @note 参考: Byrd-Lu-Nocedal-Zhu (TOMS 778)
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> lbfgsbMinimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const V& x0,
    const V& lower,
    const V& upper,
    const LBFGSBOptions<T>& options = LBFGSBOptions<T>(),
    const LineSearchParams<T>& line_search_params = LineSearchParams<T>());

// =====================================================================
// ペナルティ法（Penalty Method）
// =====================================================================

/**
 * @brief ペナルティ法による制約付き最適化
 *
 * 不等式制約 c_i(x) ≤ 0 をペナルティ項に変換し、
 * Φ(x) = f(x) + μ Σ max(0, c_i(x))² を無制約最適化する。
 * μ を段階的に増加させて制約違反を抑制。
 *
 * @param f 目的関数
 * @param df 目的関数の勾配
 * @param constraints 不等式制約関数の列 (c_i(x) ≤ 0)
 * @param x0 初期点
 * @param options 制約付き最適化オプション
 * @param line_search_params ライン探索パラメータ
 * @return 最適化結果
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> penaltyMinimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const std::vector<std::function<T(const V&)>>& constraints,
    const V& x0,
    const ConstrainedOptions<T>& options = ConstrainedOptions<T>(),
    const LineSearchParams<T>& line_search_params = LineSearchParams<T>());

// =====================================================================
// 拡張ラグランジュ法（Augmented Lagrangian Method）
// =====================================================================

/**
 * @brief 拡張ラグランジュ法による制約付き最適化
 *
 * 等式制約 h_i(x) = 0 と不等式制約 g_i(x) ≤ 0 を処理。
 * ラグランジュ乗数を推定して更新することで、
 * 純粋ペナルティ法の ill-conditioning を回避。
 *
 * L(x,λ,ν,μ) = f(x) + Σ λ_i h_i + (μ/2) Σ h_i²
 *             + (μ/2) Σ max(0, ν_i/μ + g_i)²
 *
 * @param f 目的関数
 * @param df 目的関数の勾配
 * @param eq_constraints 等式制約 h_i(x) = 0
 * @param ineq_constraints 不等式制約 g_i(x) ≤ 0
 * @param x0 初期点
 * @param options 制約付き最適化オプション
 * @param line_search_params ライン探索パラメータ
 * @return 最適化結果
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> augmentedLagrangianMinimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    const std::vector<std::function<T(const V&)>>& eq_constraints,
    const std::vector<std::function<T(const V&)>>& ineq_constraints,
    const V& x0,
    const ConstrainedOptions<T>& options = ConstrainedOptions<T>(),
    const LineSearchParams<T>& line_search_params = LineSearchParams<T>());


// =====================================================================
// Hooke-Jeeves パターン探索（微分不要の最適化）
// =====================================================================

/**
 * @brief Hooke-Jeeves パターン探索法
 *
 * 微分不要の直接探索法。各座標方向に ±dx で試行（探索移動）し、
 * 改善が見つかれば前回位置からの外挿方向にパターン移動する。
 * ステップサイズを √2 倍で拡大/縮小し、全変数の dx が min_step 以下で収束。
 *
 * @param f 最小化する関数
 * @param x0 初期点
 * @param options Hooke-Jeeves オプション
 * @return 最適化結果
 *
 * @note 出典: 今野浩・山下浩「非線形計画法」p281-282
 * @note lib++20 からの改善:
 *   - 仮想関数のクラス → std::function の自由関数
 *   - OptimizationResult 統一インターフェース
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> hookeJeevesMinimize(
    const std::function<T(const V&)>& f,
    const V& x0,
    const HookeJeevesOptions<T>& options = HookeJeevesOptions<T>());

// =====================================================================
// 遺伝的アルゴリズム（実数コード型）
// =====================================================================

/**
 * @brief 実数コード型遺伝的アルゴリズムによる大域的最適化
 *
 * 探索空間 [lower, upper] 内で関数 f の最小値を求める。
 * - 選択: トーナメント選択
 * - 交叉: BLX-α (Blend Crossover, α=0.5)
 * - 突然変異: ガウス変異
 * - エリート保存
 *
 * @param f 最小化する関数
 * @param lower 下限ベクトル
 * @param upper 上限ベクトル
 * @param options GA オプション
 * @return 最適化結果
 *
 * @note lib++20 からの改善:
 *   - ビット表現 → 実数コード（連続最適化に適する）
 *   - ルーレット選択 → トーナメント選択（選択圧の制御が容易）
 *   - 一点交叉 → BLX-α（実数変数に適する）
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> geneticAlgorithmMinimize(
    const std::function<T(const V&)>& f,
    const V& lower,
    const V& upper,
    const GeneticAlgorithmOptions<T>& options = GeneticAlgorithmOptions<T>());

// =====================================================================
// 焼きなまし法 (Simulated Annealing)
// =====================================================================

/**
 * @brief 焼きなまし法による最小化
 *
 * 確率的メタヒューリスティクス。Metropolis 判定により局所解からの脱出が可能。
 * 指数冷却スケジュール: temp(k+1) = cooling_rate * temp(k)
 *
 * @param f 最小化する関数
 * @param lower 探索範囲の下限
 * @param upper 探索範囲の上限
 * @param options SA オプション
 * @return 最適化結果
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> simulatedAnnealingMinimize(
    const std::function<T(const V&)>& f,
    const V& lower, const V& upper,
    const SimulatedAnnealingOptions<T>& options = SimulatedAnnealingOptions<T>());

// =====================================================================
// 差分進化法 (Differential Evolution, DE/rand/1/bin)
// =====================================================================

/**
 * @brief 差分進化法による最小化
 *
 * 集団ベースの確率的最適化。微分不要。
 * 戦略: DE/rand/1/bin (ランダム基底ベクトル, 1差分ベクトル, 二項交叉)
 *
 * @param f 最小化する関数
 * @param lower 探索範囲の下限
 * @param upper 探索範囲の上限
 * @param options DE オプション
 * @return 最適化結果
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> differentialEvolutionMinimize(
    const std::function<T(const V&)>& f,
    const V& lower, const V& upper,
    const DifferentialEvolutionOptions<T>& options = DifferentialEvolutionOptions<T>());

// =====================================================================
// CMA-ES (共分散行列適応進化戦略)
// =====================================================================

/**
 * @brief CMA-ES による最小化
 *
 * 確率的な進化戦略。共分散行列 C を適応的に更新して探索分布の形状を学習する。
 * 微分不要。中規模の連続最適化問題（n < 数百）に適する。
 * ill-conditioned / non-separable 問題に強い。
 *
 * 参考: Hansen & Ostermeier (2001), Hansen (2016) "The CMA Evolution Strategy: A Tutorial"
 *
 * @param f 最小化する関数
 * @param lower 探索範囲の下限
 * @param upper 探索範囲の上限
 * @param options CMA-ES オプション
 * @return 最適化結果
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> cmaesMinimize(
    const std::function<T(const V&)>& f,
    const V& lower, const V& upper,
    const CMAESOptions<T>& options = CMAESOptions<T>());

// =====================================================================
// jSO (改良版適応的差分進化法)
// =====================================================================

/**
 * @brief jSO (改良版 L-SHADE) による最小化
 *
 * SHADE の成功履歴に基づく F/CR 適応 + L-SHADE の線形人口縮小
 * + jSO 固有の F/CR 初期化改良を組み合わせた差分進化法。
 * 微分不要。大域的探索に優れる。
 *
 * 参考: Brest et al. (2017) "Single Objective Real-Parameter Optimization: Algorithm jSO"
 *
 * 主要コンポーネント:
 * - current-to-pbest/1 変異戦略
 * - 成功履歴に基づく F, CR の適応 (Lehmer 平均 / 重み付き算術平均)
 * - 線形人口サイズ縮小 (NP_init → NP_min)
 * - 外部アーカイブ (多様性維持)
 *
 * @param f 最小化する関数
 * @param lower 探索範囲の下限
 * @param upper 探索範囲の上限
 * @param options jSO オプション
 * @return 最適化結果
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> jsoMinimize(
    const std::function<T(const V&)>& f,
    const V& lower, const V& upper,
    const JSOOptions<T>& options = JSOOptions<T>());

// MKL対応の最適化アルゴリズム（MKLが利用可能な場合）
// =====================================================================
// 非線形最小二乗フィッティング (Levenberg-Marquardt)
// =====================================================================

/**
 * @brief Levenberg-Marquardt 法による非線形最小二乗フィッティング
 *
 * 残差関数 r(p) の二乗和 Σ r_i² を最小化し、
 * パラメータ共分散行列・標準誤差等の統計量を返す。
 *
 * 適応的 λ 制御 (gain ratio ベース):
 *   (J^T J + λ I) δ = -J^T r
 *   ρ = (cost - cost_new) / predicted_reduction
 *   ρ > 0 → 採択, λ *= lambdaDown
 *   ρ ≤ 0 → 棄却, λ *= lambdaUp
 *
 * @param residuals 残差関数 r(p): R^n → R^m
 * @param jacobian  ヤコビアン J(p): R^n → R^{m×n}
 * @param p0        初期パラメータ (n 次元)
 * @param options   LM オプション
 * @return フィッティング結果 (パラメータ, 共分散, 統計量)
 */
template<concepts::OrderedField T>
LeastSquaresFitResult<T> leastSquaresFit(
    const std::function<Vector<T>(const Vector<T>&)>& residuals,
    const std::function<Matrix<T>(const Vector<T>&)>& jacobian,
    const Vector<T>& p0,
    const LeastSquaresOptions<T>& options = LeastSquaresOptions<T>());

/**
 * @brief curveFit — 高レベル曲線フィッティング (scipy.optimize.curve_fit 相当)
 *
 * モデル関数 y = model(x, params) を観測データにフィットする。
 * 残差とヤコビアンを内部で自動構築し leastSquaresFit を呼ぶ。
 * ヤコビアンは中心差分で数値計算。
 *
 * @param model  モデル関数 f(x, params) → y
 * @param xdata  独立変数 (長さ m)
 * @param ydata  従属変数 (長さ m)
 * @param p0     初期パラメータ (n 次元)
 * @param options LM オプション
 * @return フィッティング結果
 */
template<concepts::OrderedField T>
LeastSquaresFitResult<T> curveFit(
    const std::function<T(T, const Vector<T>&)>& model,
    const std::vector<T>& xdata,
    const std::vector<T>& ydata,
    const Vector<T>& p0,
    const LeastSquaresOptions<T>& options = LeastSquaresOptions<T>());

} // namespace sangi

// 実装は optimization_nd_impl.hpp に分離。
// src/math/optimization/optimization.cpp で float/double の明示的インスタンス化を行う。

#endif // SANGI_OPTIMIZATION_ND_HPP