// 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 Backtracking line search
 * @tparam T Ordered field type (such as the reals)
 * @tparam V Vector type
 * @param f Function to minimize
 * @param x Current point
 * @param direction Search direction
 * @param fx Function value at the current point
 * @param gradient Gradient at the current point
 * @param params Line search parameters
 * @return Tuple of (alpha, new point, new function value, success flag)
 */
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 Line search (steepest descent method)
 * @tparam T Ordered field type (such as the reals)
 * @tparam V Vector type
 * @param f Function to minimize
 * @param df Function that computes the gradient of the function
 * @param x0 Initial point
 * @param options Optimization options
 * @param line_search_params Line search parameters
 * @return Tuple of the minimizing point, the minimum value, and a success flag
 */
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 Quasi-Newton method (BFGS method)
 * @tparam T Ordered field type (such as the reals)
 * @tparam V Vector type
 * @tparam M Matrix type
 * @param f Function to minimize
 * @param df Function that computes the gradient of the function
 * @param x0 Initial point
 * @param options Optimization options
 * @param line_search_params Line search parameters
 * @return Tuple of the minimizing point, the minimum value, and a success flag
 */
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 Newton conjugate gradient method (Truncated Newton Method)
 * @tparam T Ordered field type (such as the reals)
 * @tparam V Vector type
 * @param f Function to minimize
 * @param df Function that computes the gradient of the function
 * @param d2f Function that computes the product of the function's Hessian matrix and a direction vector
 * @param x0 Initial point
 * @param options Optimization options
 * @param line_search_params Line search parameters
 * @return Tuple of the minimizing point, the minimum value, and a success flag
 */
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 Outer product of vectors (the result is a matrix)
 * @tparam T Field type (such as the reals)
 * @tparam V Vector type
 * @tparam M Matrix type
 * @param a Vector
 * @param b Vector
 * @return Outer product a⊗b (matrix with entries 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 simplex method (derivative-free optimization)
// =====================================================================

/**
 * @brief Nelder-Mead simplex method
 * @tparam T Ordered field type
 * @tparam V Vector type
 * @param f Function to minimize
 * @param x0 Initial point
 * @param options Nelder-Mead options
 * @return Optimization result
 *
 * @note Improvements over sangi:
 *   - Convergence test: dual test on function-value spread AND simplex diameter
 *   - On timeout: returns the best point (sangi returned TYPE(0))
 *   - Standard Nelder-Mead coefficients (α=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 method (nonlinear least squares)
// =====================================================================

/**
 * @brief Gauss-Newton method (with Levenberg-Marquardt regularization)
 * @tparam T Ordered field type
 * @tparam V Vector type
 * @tparam M Matrix type
 * @param residuals Residual function r(x): V→V
 * @param jacobian Jacobian matrix function J(x): V→M
 * @param x0 Initial parameters
 * @param options Gauss-Newton options
 * @param line_search_params Line search parameters
 * @return Optimization result (value is the sum of squared residuals)
 *
 * @note Improvements over sangi:
 *   - LM damping (J^T J + λI) avoids singularity
 *   - Removed fallback to the old dX → reported as a failure instead
 *   - Added Armijo line search (no longer a fixed step)
 *   - Removed stdout debug output
 */
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>());

// =====================================================================
// Nonlinear conjugate gradient method (Fletcher-Reeves / Polak-Ribière)
// =====================================================================

/**
 * @brief Nonlinear conjugate gradient method
 * @tparam T Ordered field type
 * @tparam V Vector type
 * @param f Function to minimize
 * @param df Gradient function
 * @param x0 Initial point
 * @param options Conjugate gradient options
 * @param line_search_params Line search parameters
 * @return Optimization result
 *
 * @note Improvements over sangi:
 *   - Armijo backtracking (sangi used a fixed α=0.000001 + unbounded bracket expansion)
 *   - Restart every n steps (sangi left loss of conjugacy unaddressed)
 *   - Zero-division guard in the β computation
 */
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 method (memory-efficient quasi-Newton method)
 *
 * A memory-efficient version of BFGS. Instead of explicitly storing the inverse
 * Hessian, it computes H⁻¹g in O(mn) via the two-loop recursion using the past
 * m gradient differences. Standard for large-scale problems (1000+ variables).
 *
 * @param f Function to minimize
 * @param df Gradient function
 * @param x0 Initial point
 * @param options L-BFGS options (memory_size = number of vectors to retain)
 * @param line_search_params Line search parameters
 * @return Optimization result
 */
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 (bound-constrained L-BFGS)
// =====================================================================

/**
 * @brief L-BFGS-B method (bound-constrained memory-efficient quasi-Newton method)
 *
 * Optimization where each variable has upper and lower bound constraints l_i ≤ x_i ≤ u_i.
 * It projects the search direction onto the feasible region via the gradient
 * projection method, and computes the L-BFGS step using only the free variables.
 *
 * @param f Function to minimize
 * @param df Gradient function
 * @param x0 Initial point (automatically projected if not feasible)
 * @param lower Lower bound (-numeric_limits::max() for -inf)
 * @param upper Upper bound (+numeric_limits::max() for +inf)
 * @param options L-BFGS-B options
 * @param line_search_params Line search parameters
 * @return Optimization result
 *
 * @note Reference: 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 Constrained optimization via the penalty method
 *
 * Converts inequality constraints c_i(x) ≤ 0 into penalty terms and performs
 * unconstrained optimization of Φ(x) = f(x) + μ Σ max(0, c_i(x))².
 * μ is increased gradually to suppress constraint violations.
 *
 * @param f Objective function
 * @param df Gradient of the objective function
 * @param constraints Sequence of inequality constraint functions (c_i(x) ≤ 0)
 * @param x0 Initial point
 * @param options Constrained optimization options
 * @param line_search_params Line search parameters
 * @return Optimization result
 */
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 Constrained optimization via the augmented Lagrangian method
 *
 * Handles equality constraints h_i(x) = 0 and inequality constraints g_i(x) ≤ 0.
 * By estimating and updating the Lagrange multipliers, it avoids the
 * ill-conditioning of the pure penalty method.
 *
 * L(x,λ,ν,μ) = f(x) + Σ λ_i h_i + (μ/2) Σ h_i²
 *             + (μ/2) Σ max(0, ν_i/μ + g_i)²
 *
 * @param f Objective function
 * @param df Gradient of the objective function
 * @param eq_constraints Equality constraints h_i(x) = 0
 * @param ineq_constraints Inequality constraints g_i(x) ≤ 0
 * @param x0 Initial point
 * @param options Constrained optimization options
 * @param line_search_params Line search parameters
 * @return Optimization result
 */
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>());

// =====================================================================
// Sequential quadratic programming (SQP)
// =====================================================================

namespace detail {

// Small-scale QP solver (active set method)
// min 0.5 x^T H x + c^T x  s.t. A_eq x = b_eq, A_ineq x ≤ b_ineq
template<concepts::OrderedField T>
Vector<T> solveQP(
    const Matrix<T>& H,
    const Vector<T>& c,
    const Matrix<T>& A_eq,
    const Vector<T>& b_eq,
    const Matrix<T>& A_ineq,
    const Vector<T>& b_ineq,
    const Vector<T>& x0,
    std::size_t max_iter = 200);

} // namespace detail

/**
 * @brief Constrained optimization via sequential quadratic programming (SQP)
 *
 * At each iteration it makes a quadratic approximation of the objective function
 * and a linear approximation of the constraints, then solves the QP subproblem.
 * The standard method for nonlinear constrained optimization.
 * Uses BFGS updates for the Hessian approximation.
 *
 * @param f Objective function
 * @param df Gradient of the objective function
 * @param eq_constraints Equality constraints h_i(x) = 0
 * @param ineq_constraints Inequality constraints g_i(x) ≤ 0
 * @param x0 Initial point
 * @param options Constrained optimization options
 * @return Optimization result
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> sqpMinimize(
    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>());

// =====================================================================
// Interior Point Method
// =====================================================================

/**
 * @brief Inequality-constrained optimization via the interior point method (barrier method)
 *
 * Handles inequality constraints g_i(x) ≤ 0 with a barrier function:
 * B(x,t) = f(x) - t Σ log(-g_i(x))
 * As t is decreased toward 0, it approaches the constraint boundary from the inside.
 *
 * @param f Objective function
 * @param df Gradient of the objective function
 * @param ineq_constraints Inequality constraints g_i(x) ≤ 0
 * @param x0 Initial point (an interior point strictly satisfying all constraints, g_i(x0) < 0)
 * @param options Constrained optimization options
 * @param line_search_params Line search parameters
 * @return Optimization result
 * @throw std::invalid_argument If the initial point is not feasible
 */
template<concepts::OrderedField T, typename V>
    requires concepts::VectorOf<V, T>
[[nodiscard]] OptimizationResult<V, T> interiorPointMinimize(
    const std::function<T(const V&)>& f,
    const std::function<V(const V&)>& df,
    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 pattern search (derivative-free optimization)
// =====================================================================

/**
 * @brief Hooke-Jeeves pattern search method
 *
 * A derivative-free direct search method. It probes ±dx along each coordinate
 * direction (exploratory move), and if an improvement is found, makes a pattern
 * move in the extrapolated direction from the previous position.
 * The step size is expanded/contracted by a factor of √2, and convergence is
 * reached when dx for all variables falls below min_step.
 *
 * @param f Function to minimize
 * @param x0 Initial point
 * @param options Hooke-Jeeves options
 * @return Optimization result
 *
 * @note Source: Hiroshi Konno & Hiroshi Yamashita, "Nonlinear Programming," pp. 281-282
 * @note Improvements over sangi:
 *   - Virtual-function class → free function using std::function
 *   - Unified OptimizationResult interface
 */
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>());

// =====================================================================
// Genetic algorithm (real-coded)
// =====================================================================

/**
 * @brief Global optimization via a real-coded genetic algorithm
 *
 * Finds the minimum of function f within the search space [lower, upper].
 * - Selection: tournament selection
 * - Crossover: BLX-α (Blend Crossover, α=0.5)
 * - Mutation: Gaussian mutation
 * - Elitism
 *
 * @param f Function to minimize
 * @param lower Lower bound vector
 * @param upper Upper bound vector
 * @param options GA options
 * @return Optimization result
 *
 * @note Improvements over sangi:
 *   - Bit representation → real coding (suited to continuous optimization)
 *   - Roulette selection → tournament selection (easier control of selection pressure)
 *   - Single-point crossover → BLX-α (suited to real-valued variables)
 */
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 Minimization via simulated annealing
 *
 * A stochastic metaheuristic. The Metropolis criterion allows escaping from local optima.
 * Exponential cooling schedule: temp(k+1) = cooling_rate * temp(k)
 *
 * @param f Function to minimize
 * @param lower Lower bound of the search range
 * @param upper Upper bound of the search range
 * @param options SA options
 * @return Optimization result
 */
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 Minimization via differential evolution
 *
 * Population-based stochastic optimization. Derivative-free.
 * Strategy: DE/rand/1/bin (random base vector, 1 difference vector, binomial crossover)
 *
 * @param f Function to minimize
 * @param lower Lower bound of the search range
 * @param upper Upper bound of the search range
 * @param options DE options
 * @return Optimization result
 */
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 (Covariance Matrix Adaptation Evolution Strategy)
// =====================================================================

/**
 * @brief Minimization via CMA-ES
 *
 * A stochastic evolution strategy. It adaptively updates the covariance matrix C
 * to learn the shape of the search distribution.
 * Derivative-free. Suited to medium-scale continuous optimization problems (n < a few hundred).
 * Strong on ill-conditioned / non-separable problems.
 *
 * Reference: Hansen & Ostermeier (2001), Hansen (2016) "The CMA Evolution Strategy: A Tutorial"
 *
 * @param f Function to minimize
 * @param lower Lower bound of the search range
 * @param upper Upper bound of the search range
 * @param options CMA-ES options
 * @return Optimization result
 */
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 (improved adaptive differential evolution)
// =====================================================================

/**
 * @brief Minimization via jSO (improved L-SHADE)
 *
 * A differential evolution method combining SHADE's success-history-based F/CR
 * adaptation + L-SHADE's linear population reduction
 * + jSO's own improved F/CR initialization.
 * Derivative-free. Excellent at global search.
 *
 * Reference: Brest et al. (2017) "Single Objective Real-Parameter Optimization: Algorithm jSO"
 *
 * Key components:
 * - current-to-pbest/1 mutation strategy
 * - Success-history-based adaptation of F, CR (Lehmer mean / weighted arithmetic mean)
 * - Linear population size reduction (NP_init → NP_min)
 * - External archive (diversity maintenance)
 *
 * @param f Function to minimize
 * @param lower Lower bound of the search range
 * @param upper Upper bound of the search range
 * @param options jSO options
 * @return Optimization result
 */
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-enabled optimization algorithms (when MKL is available)
// =====================================================================
// Nonlinear least-squares fitting (Levenberg-Marquardt)
// =====================================================================

/**
 * @brief Nonlinear least-squares fitting via the Levenberg-Marquardt method
 *
 * Minimizes the sum of squares Σ r_i² of the residual function r(p), and returns
 * statistics such as the parameter covariance matrix and standard errors.
 *
 * Adaptive λ control (gain-ratio based):
 *   (J^T J + λ I) δ = -J^T r
 *   ρ = (cost - cost_new) / predicted_reduction
 *   ρ > 0 → accept, λ *= lambdaDown
 *   ρ ≤ 0 → reject, λ *= lambdaUp
 *
 * @param residuals Residual function r(p): R^n → R^m
 * @param jacobian  Jacobian J(p): R^n → R^{m×n}
 * @param p0        Initial parameters (n-dimensional)
 * @param options   LM options
 * @return Fitting result (parameters, covariance, statistics)
 */
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 — high-level curve fitting (equivalent to scipy.optimize.curve_fit)
 *
 * Fits the model function y = model(x, params) to observed data.
 * Internally builds the residuals and Jacobian automatically and calls leastSquaresFit.
 * The Jacobian is computed numerically by central differences.
 *
 * @param model  Model function f(x, params) → y
 * @param xdata  Independent variable (length m)
 * @param ydata  Dependent variable (length m)
 * @param p0     Initial parameters (n-dimensional)
 * @param options LM options
 * @return Fitting result
 */
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

// The implementation is separated into optimization_nd_impl.hpp.
// Explicit instantiation for float/double is done in src/math/optimization/optimization.cpp.

#endif // SANGI_OPTIMIZATION_ND_HPP