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

// SimulatedAnnealing.hpp — Simulated Annealing
//
// Overview:
//   Simulated annealing (SA) is a metaheuristic applicable to both combinatorial
//   and continuous optimization. Its hallmark is the ability to escape local
//   optima probabilistically.
//
// Algorithm outline:
//   1. Generate a neighbor x' from the current solution x (neighbor function)
//   2. Compute ΔE = E(x') - E(x)
//   3. Accept if ΔE < 0 (improvement). Even when ΔE >= 0, accept with probability exp(-ΔE/T) (Metropolis criterion)
//   4. Cool the temperature T by coolingRate (exponential cooling schedule)
//   5. Terminate when T < tempMin or maxIter is reached
//
// Tuning tips:
//   - tempInit:    set to roughly 5-10 times the typical variation of the objective function
//   - coolingRate: 0.99 to 0.9999 (trade-off between problem size and accuracy)
//   - maxIter:     take it large enough. With coolingRate=0.995, reaching tempInit/tempMin
//                  requires about -log(tempMin/tempInit)/log(1/coolingRate) iterations
//
// Usage example:
//   // 1D: minimize f(x) = (x-3)^2 (over the range [-10, 10])
//   auto energy = [](const std::vector<double>& x) { return (x[0]-3)*(x[0]-3); };
//   auto neighbor = [](const std::vector<double>& x, std::mt19937& rng) {
//       std::uniform_real_distribution<double> d(-0.5, 0.5);
//       return std::vector<double>{ x[0] + d(rng) };
//   };
//   auto result = sangi::simulatedAnnealing<double>(energy, neighbor, {0.0});
//   // result.bestSolution[0] ≈ 3.0
//   // result.bestEnergy      ≈ 0.0
//   // result.iterations      = actual number of iterations

#ifndef SANGI_OPTIMIZATION_SIMULATED_ANNEALING_HPP
#define SANGI_OPTIMIZATION_SIMULATED_ANNEALING_HPP

#include <vector>
#include <functional>
#include <random>
#include <cmath>

namespace sangi {

    // Result of simulated annealing
    template<typename T>
    struct SAResult {
        std::vector<T> bestSolution;  // best solution found over the entire history
        T bestEnergy;                  // energy value of the best solution
        size_t iterations;             // actual number of iterations (≤ maxIter)
    };

    // simulatedAnnealing — minimization by simulated annealing
    //
    //   energy:      objective function E(x): std::vector<T> → T  (to be minimized)
    //   neighbor:    neighbor generation function (x, rng) → x'  (determines the search step size)
    //   x0:          initial solution vector
    //   tempInit:    initial temperature (larger means broader initial search)
    //   tempMin:     minimum temperature (terminate once below this)
    //   coolingRate: cooling rate 0 < coolingRate < 1  (closer to 1 means slower cooling)
    //   maxIter:     maximum number of iterations (upper bound even if tempMin is not reached)
    //   seed:        random seed (can be fixed for reproducibility)
    template<typename T>
    SAResult<T> simulatedAnnealing(
        std::function<T(const std::vector<T>&)> energy,
        std::function<std::vector<T>(const std::vector<T>&, std::mt19937&)> neighbor,
        const std::vector<T>& x0,
        T tempInit = T(100), T tempMin = T(1e-8), T coolingRate = T(0.995),
        size_t maxIter = 100000, unsigned seed = 42)
    {
        std::mt19937 rng(seed);
        std::uniform_real_distribution<T> uniform(T(0), T(1));

        std::vector<T> current = x0;
        T currentEnergy = energy(current);

        std::vector<T> best = current;
        T bestEnergy = currentEnergy;

        T temp = tempInit;
        size_t iter = 0;

        for (; iter < maxIter && temp > tempMin; ++iter) {
            std::vector<T> candidate = neighbor(current, rng);
            T candidateEnergy = energy(candidate);
            T delta = candidateEnergy - currentEnergy;

            // Metropolis test: accept if energy decreases, or with probability exp(-ΔE/T)
            if (delta < T(0) || uniform(rng) < std::exp(-delta / temp)) {
                current = std::move(candidate);
                currentEnergy = candidateEnergy;

                if (currentEnergy < bestEnergy) {
                    best = current;
                    bestEnergy = currentEnergy;
                }
            }

            temp *= coolingRate;
        }

        return SAResult<T>{ std::move(best), bestEnergy, iter };
    }

} // namespace sangi

#endif // SANGI_OPTIMIZATION_SIMULATED_ANNEALING_HPP