genetic_algorithm.h

/* A C++20 implementation of the Genetic Algorithm

   Based on scikit-opt's GA implementation (guofei9987/scikit-opt)
   Rewritten in modern C++20 for the ExOpt framework.

   Design: GaSolver does NOT own the objective function. It generates
   candidate positions, accepts externally-computed fitness values, and
   advances the population by one generation per Update() call.

   The solver supports two encoding modes:
     - Binary:  Gray-coded binary strings (classic GA)
     - Real:    real-valued chromosomes (RCGA-style)

   Licensed under the MIT License.
*/

#ifndef EX_OPT_GA_H_
#define EX_OPT_GA_H_

#include <algorithm>
#include <cmath>
#include <functional>
#include <iostream>
#include <limits>
#include <numeric>
#include <random>
#include <span>
#include <vector>

#include "policy/base_policy.h"

namespace ex_opt::ga {

// ============================================================
// CONSTANTS
// ============================================================
inline constexpr int kDefaultPopSize = 50;
inline constexpr int kDefaultMaxIter = 200;
inline constexpr double kDefaultProbMut = 0.001;
inline constexpr double kDefaultProbCros = 0.9;

// ============================================================
// STRATEGY ENUMS
// ============================================================
enum class EncodingType {
  Binary,  // Gray-coded binary strings → continuous variables
  Real     // Real-valued chromosomes (RCGA)
};

enum class SelectionType {
  Tournament,  // fast vectorized tournament selection
  Roulette     // fitness-proportional roulette wheel
};

enum class CrossoverType {
  TwoPoint,  // two-point crossover
  OnePoint,  // single-point crossover
  Pmx        // PMX for permutation encoding
};

// ============================================================
// SETTINGS
// ============================================================
struct GaSettings {
  SearchSpace space;  // search-space definition (satisfies Boundary)

  double goal = 1e-7;  // error threshold for early stop

  int size_pop = kDefaultPopSize;  // population size (must be even)
  int max_iter = kDefaultMaxIter;  // maximum generations
  int early_stop = 0;              // consecutive gens w/o improvement before stop (0 = off)
  int print_every = 100;           // logging interval (0 = silent)
  int step = 0;                    // current generation (updated by solver)

  // ---- Encoding ----
  EncodingType encoding = EncodingType::Binary;
  double precision = 1e-7;  // per-variable precision (binary encoding only)

  // ---- Operator probabilities ----
  double prob_mut = kDefaultProbMut;    // per-gene mutation probability
  double prob_cros = kDefaultProbCros;  // crossover probability (RCGA only)

  // ---- Selection ----
  SelectionType selection = SelectionType::Tournament;
  int tourn_size = 3;  // tournament size

  // ---- Crossover ----
  CrossoverType crossover = CrossoverType::TwoPoint;

  /// Construct with uniform bounds for all dimensions.
  GaSettings(int dim, double range_lo, double range_hi);

  // ---- Boundary accessors ----
  [[nodiscard]] int dim() const { return space.dim(); }
  [[nodiscard]] double lower(int i) const { return space.lower(i); }
  [[nodiscard]] double upper(int i) const { return space.upper(i); }
};

static_assert(Boundary<GaSettings>);

// ============================================================
// RESULT
// ============================================================
struct GaResult {
  double error = std::numeric_limits<double>::max();
  std::vector<double> best_x;
  std::vector<double> generation_best_y;  // best Y per generation
};

// ============================================================
// GA SOLVER
//
//   Lifecycle:
//     1. Initialize()          — create random initial population
//     2. GetPositions()        — expose decoded positions for
//                                 external evaluation
//     3. Update(fitness)       — feed fitness back; solver runs
//                                 selection → crossover → mutation,
//                                 then decodes the next population
//     4. repeat from 2 until Done()
// ============================================================
class GaSolver {
 public:
  explicit GaSolver(GaSettings settings);

  /// Create random initial population and decode to X.
  /// Must be called once before Update().
  void Initialize();

  /// All decoded candidate positions as a flat span
  /// [individual * dim + d]. Length = GetPopSize() * GetDimensions().
  [[nodiscard]] std::span<const double> GetPositions() const { return X_; }

  /// Submit fitness values for the *current* positions (one per
  /// individual), then run one full GA generation:
  ///   ranking → selection → crossover → mutation → decode
  ///
  /// After this call, GetPositions() returns the *next* generation
  /// of decoded positions, ready for another round of evaluation.
  void Update(std::span<const double> fitness);

  /// True when goal reached, max generations exhausted, or
  /// early_stop triggered.
  [[nodiscard]] bool Done() const { return done_; }

  [[nodiscard]] int GetCurrentGeneration() const { return current_gen_; }
  [[nodiscard]] int GetMaxGenerations() const { return settings_.max_iter; }
  [[nodiscard]] int GetPopSize() const { return settings_.size_pop; }
  [[nodiscard]] int GetDimensions() const { return settings_.dim(); }
  [[nodiscard]] const GaResult& GetResult() const { return result_; }
  [[nodiscard]] const GaSettings& GetSettings() const { return settings_; }

 private:
  GaSettings settings_;
  std::mt19937 rng_;

  // ---- Population state ----
  std::vector<int> Chrom_;         // binary chromosomes  [pop * len_chrom_]
  std::vector<double> ChromReal_;  // real chromosomes   [pop * dim]
  std::vector<double> X_;          // decoded positions  [pop * dim]

  // ---- Per-generation working state ----
  std::vector<double> Y_;     // fitness (lower = better)
  std::vector<double> FitV_;  // scaled fitness (higher = better)

  // ---- Binary-encoding metadata ----
  std::vector<int> Lind_;  // bits per variable
  int len_chrom_ = 0;      // total chromosome length in bits

  GaResult result_;
  int current_gen_ = 0;
  bool initialized_ = false;
  bool done_ = false;

  // ---- Initialisation helpers ----
  void crtbp_binary_();
  void crtbp_real_();
  void chrom2x_binary_();
  void chrom2x_real_();

  /// Gray code → [0, 1) real value for a single variable segment.
  /// `base` points to the first bit of the segment in Chrom_.
  [[nodiscard]] double gray2rv_(int base, int lind) const;

  // ---- Ranking ----
  void ranking_();

  // ---- Selection ----
  void selection_tournament_();
  void selection_roulette_();

  // ---- Crossover ----
  void crossover_2point_binary_();
  void crossover_1point_binary_();
  void crossover_sbx_();  // simulated binary crossover (real encoding)

  // ---- Mutation ----
  void mutate_bitflip_();
  void mutate_polynomial_();

  // ---- Random utilities ----
  [[nodiscard]] double uniform_();
  [[nodiscard]] int uniform_int_(int max);
};

// ============================================================
// GENETIC ALGORITHM POLICY  (implements BatchPolicy<State>)
// ============================================================
class GeneticAlgorithmPolicy : public BatchPolicy<std::vector<double>> {
 public:
  explicit GeneticAlgorithmPolicy(GaSettings settings) : solver_(std::move(settings)) {}

  /// Create the initial random population.
  void Initialize() override { solver_.Initialize(); }

  /// Number of candidate positions to evaluate per generation.
  [[nodiscard]] int BatchSize() const override { return solver_.GetPopSize(); }

  /// Extract the i-th candidate as a std::vector<double>.
  [[nodiscard]] std::vector<double> GetCandidate(int i) const override {
    const int d = solver_.GetDimensions();
    auto pos = solver_.GetPositions().subspan(i * d, d);
    std::vector<double> s(d);
    std::ranges::copy(pos, s.begin());
    return s;
  }

  /// Feed fitness values (one per individual) back to the GA.
  /// The solver runs one generation and decodes the next population.
  void SubmitFitness(std::span<const double> fitness) override {
    solver_.Update(fitness);
    // Sync best state
    best_state_.resize(solver_.GetDimensions());
    std::ranges::copy(solver_.GetResult().best_x, best_state_.begin());
    done_ = solver_.Done();
  }

  [[nodiscard]] std::vector<double> GetBestState() const override { return best_state_; }
  [[nodiscard]] bool Done() const override { return done_; }

  /// Expose the underlying solver for fine-grained control.
  [[nodiscard]] auto& GetSolver() { return solver_; }
  [[nodiscard]] const auto& GetSolver() const { return solver_; }

 private:
  GaSolver solver_;
  std::vector<double> best_state_;
  bool done_ = false;
};

}  // namespace ex_opt::ga

#endif  // EX_OPT_GA_H_