// 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 {

// Common options struct shared by optimization algorithms
template<concepts::OrderedField T>
struct OptimizationOptions {
    T tolerance = std::numeric_limits<T>::epsilon() * 100;
    size_t max_iterations = 1000;
    bool verbose = false;
};

// Optimization result struct
template<typename PointType, concepts::OrderedField ValueType>
struct OptimizationResult {
    PointType point;       // optimized point
    ValueType value;       // optimized value
    bool success;          // optimization success flag
    size_t iterations;     // number of iterations performed

    // Static helper to create a success result
    static OptimizationResult makeSuccess(const PointType& p, const ValueType& v, size_t iters) {
        return {p, v, true, iters};
    }

    // Static helper to create a failure result
    static OptimizationResult makeFailure(const PointType& p, const ValueType& v, size_t iters) {
        return {p, v, false, iters};
    }
};

// Result type for one-dimensional problems
template<concepts::OrderedField T>
using OptimizationResult1D = OptimizationResult<T, T>;

// Line-search parameters struct
template<concepts::OrderedField T>
struct LineSearchParams {
    T alpha_init = T(0.1);          // initial step size
    T armijo_const = T(1e-4);       // constant for the Armijo rule
    T reduction_factor = T(0.5);    // backtracking reduction factor
    size_t max_iterations = 20;     // maximum number of iterations
};

// Options for the Nelder-Mead method
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);          // step size of the initial simplex
    T reflection = T(1);            // reflection coefficient α
    T expansion = T(2);             // expansion coefficient γ
    T contraction = T(0.5);         // contraction coefficient ρ
    T shrink = T(0.5);              // overall shrink coefficient σ
};

// Options for the Gauss-Newton method
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 regularization parameter λ
    bool use_line_search = true;    // whether to use Armijo line search
};

// Options for the conjugate gradient method
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 = restart every n (the dimension)
    bool use_polak_ribiere = true;  // true: Polak-Ribière, false: Fletcher-Reeves
};

// Options for 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;         // number m of past vectors to retain
};

// Options for L-BFGS-B (with bound constraints)
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;
};

// Options for Hooke-Jeeves pattern search
template<concepts::OrderedField T>
struct HookeJeevesOptions {
    T tolerance = T(1e-12);          // convergence test on change of objective value
    size_t max_iterations = 10000;   // maximum number of evaluations
    T initial_step = T(1);           // initial step size
    T min_step = T(1e-10);           // converge when the step falls to this or below
};

// Common options for constrained optimization
template<concepts::OrderedField T>
struct ConstrainedOptions {
    T tolerance = std::numeric_limits<T>::epsilon() * T(100);
    size_t max_iterations = 100;         // outer loop
    size_t inner_max_iterations = 1000;  // inner solver
    T initial_penalty = T(1);            // initial penalty parameter μ
    T penalty_growth = T(10);            // growth rate of μ
};

// Options for the genetic algorithm (real-coded)
template<concepts::OrderedField T>
struct GeneticAlgorithmOptions {
    T tolerance = T(1e-8);           // converge when change of best individual is at or below tolerance
    size_t max_generations = 1000;   // maximum number of generations
    size_t population_size = 50;     // population size
    T crossover_rate = T(0.8);       // crossover rate
    T mutation_rate = T(0.1);        // mutation rate
    T mutation_scale = T(0.1);       // mutation standard deviation (ratio relative to the search range)
    size_t tournament_size = 3;      // tournament selection size
    size_t elite_count = 2;          // number of elites preserved
    unsigned int seed = 42;          // random seed
};

// Options for simulated annealing
template<concepts::OrderedField T>
struct SimulatedAnnealingOptions {
    T tolerance = T(1e-10);              // convergence test on change of best value
    size_t max_iterations = 100000;      // maximum number of iterations
    T initial_temperature = T(100);      // initial temperature
    T cooling_rate = T(0.995);           // cooling rate: T(k+1) = cooling_rate * T(k)
    T neighbor_scale = T(0.1);           // standard deviation of Gaussian noise for neighbor generation (range ratio)
    unsigned int seed = 42;              // random seed
};

// Options for differential evolution
template<concepts::OrderedField T>
struct DifferentialEvolutionOptions {
    T tolerance = T(1e-8);               // converge when change of best individual is at or below tolerance
    size_t max_generations = 1000;       // maximum number of generations
    size_t population_size = 50;         // population size
    T mutation_factor = T(0.8);          // differential weight F
    T crossover_rate = T(0.9);           // crossover probability CR
    unsigned int seed = 42;              // random seed
};

// Options for nonlinear least-squares fitting (Levenberg-Marquardt)
template<concepts::OrderedField T>
struct LeastSquaresOptions {
    T ftol = T(1e-10);              // relative change of the residual sum of squares
    T xtol = T(1e-10);             // relative change of the parameters
    T gtol = T(1e-10);             // convergence test on the gradient norm
    std::size_t maxIterations = 200;
    T lambdaInit = T(1e-3);        // initial LM parameter λ
    T lambdaUp = T(10);            // λ growth factor (when the step is rejected)
    T lambdaDown = T(0.1);         // λ reduction factor (when the step is accepted)
};

// Result of nonlinear least-squares fitting (with statistics)
template<concepts::OrderedField T>
struct LeastSquaresFitResult {
    Vector<T> parameters;          // optimal parameters
    Vector<T> residuals;           // final residual vector
    Matrix<T> covariance;          // parameter covariance matrix s² (J^T J)^{-1}
    Vector<T> standardErrors;      // standard error of each parameter
    T chiSquared;                  // residual sum of squares Σ r_i²
    T reducedChiSquared;           // χ² / (m - n)
    std::size_t iterations;        // number of iterations
    bool success;                  // convergence flag
    std::size_t dataPoints;        // number of data points m
    std::size_t numParameters;     // number of parameters n
};

// Options for CMA-ES (Covariance Matrix Adaptation Evolution Strategy)
template<concepts::OrderedField T>
struct CMAESOptions {
    T tolerance = T(1e-10);              // convergence test on change of function value (TolFun)
    T tolx = T(1e-12);                   // convergence test on change of the mean vector (TolX)
    size_t max_generations = 10000;      // maximum number of generations
    size_t population_size = 0;          // λ (0 = automatic: 4 + floor(3*ln(n)))
    T initial_sigma = T(0.3);            // initial step size σ (ratio relative to the search range)
    T condition_limit = T(1e14);         // upper bound on the condition number of the covariance matrix (ConditionCov)
    unsigned int seed = 42;              // random seed
};

// Options for jSO (an improved adaptive differential evolution method)
template<concepts::OrderedField T>
struct JSOOptions {
    T tolerance = T(1e-8);               // converge when change of best individual is at or below tolerance
    size_t max_evaluations = 0;          // maximum number of function evaluations (0 = automatic: 10000*D)
    size_t population_size = 0;          // NP_init (0 = automatic: 25*D)
    size_t min_population_size = 4;      // NP_min (minimum population size)
    size_t history_size = 5;             // length H of the successful-parameter memory
    T archive_rate = T(2.6);             // archive size multiplier (NP_init * archive_rate)
    T p_best_rate = T(0.11);             // top fraction for pbest selection
    unsigned int seed = 42;              // random seed
};

// Check for MKL availability
inline constexpr bool has_mkl_support = SANGI_HAS_MKL == 1;

} // namespace sangi

#endif // SANGI_OPTIMIZATION_BASE_HPP