// ch13_genetic_algorithm.cpp
// Companion code for Chapter 13: Optimization and Data Fitting.
// Standalone C++11; standard library only; console output without plotting.
// Compile: c++ -std=c++11 -O2 -Wall -Wextra -pedantic ch13_genetic_algorithm.cpp -o ch13_genetic_algorithm
// Edit the example settings in main(); routines above main() can be reused.
// Fixed seeds reproduce runs with the same standard-library implementation.
// Random distributions need not give identical sequences across compilers/languages.

#include <algorithm>
#include <array>
#include <cmath>
#include <functional>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <random>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

using Vec = std::vector<double>;
using Mat = std::vector<Vec>;
using Function = std::function<double(const Vec&)>;
using ScalarFunction = std::function<double(double)>;
const double pi = std::acos(-1.0);

void require(bool condition, const std::string& message) {
    if (!condition) throw std::invalid_argument(message);
}
void finite_vector(const Vec& x) {
    require(!x.empty(), "A nonempty vector is required.");
    for (double v : x) require(std::isfinite(v), "Data must be finite.");
}
void check_bounds(const Vec& lo, const Vec& hi) {
    finite_vector(lo); finite_vector(hi);
    require(lo.size() == hi.size(), "Bounds must have matching dimensions.");
    for (std::size_t j=0; j<lo.size(); ++j)
        require(hi[j]>lo[j] && std::isfinite(hi[j]-lo[j]), "Require finite upper > lower.");
}
double objective_value(double value) {
    require(std::isfinite(value), "Objective must return a finite real scalar.");
    return value;
}
double clip(double x, double lo, double hi) { return std::max(lo,std::min(hi,x)); }
Vec linspace(double a, double b, int n) {
    require(n>=2, "Grid needs at least two points.");
    Vec x(n);
    for (int i=0; i<n; ++i) x[i]=a+(b-a)*i/(n-1);
    return x;
}
struct Random {
    std::mt19937 engine;
    std::uniform_real_distribution<double> uniform{0.0,1.0};
    std::normal_distribution<double> normal{0.0,1.0};
    explicit Random(unsigned int seed) : engine(seed) {}
    double u() { return uniform(engine); }
    double z() { return normal(engine); }
    int index(int n) { return std::uniform_int_distribution<int>(0,n-1)(engine); }
};
struct Result {
    Vec point;
    double value=0.0;
    long long evaluations=0;
    int iterations=0;
    std::string status;
};
void print_result(const Result& r) {
    std::cout << std::setprecision(12);
    for (std::size_t j=0; j<r.point.size(); ++j)
        std::cout << "p(" << j+1 << ") = " << r.point[j] << '\n';
    std::cout << "Objective value: " << r.value << "\nIterations/updates: " << r.iterations
              << "\nFunction calls: " << r.evaluations << "\nStatus: " << r.status << '\n';
}

// Real-valued chromosomes, tournament selection, uniform crossover, Gaussian
// mutation of whole children, clipping, and one unchanged elite per generation.
Result genetic_search(const Function& f, const Vec& lo, const Vec& hi, Random& rng,
                      int n, int generations, double pc, double pm, const Vec& sigma, int nt) {
    check_bounds(lo,hi); finite_vector(sigma); std::size_t d=lo.size();
    require(n>=2 && generations>0 && nt>0 && sigma.size()==d && pc>=0 && pc<=1
            && pm>=0 && pm<=1, "Invalid genetic settings.");
    for (double s:sigma) require(s>0, "Mutation widths must be positive.");
    Mat pop(n,Vec(d)); Vec values(n);
    for (int i=0;i<n;++i) {
        for (std::size_t j=0;j<d;++j) pop[i][j]=lo[j]+(hi[j]-lo[j])*rng.u();
        values[i]=objective_value(f(pop[i]));
    }
    auto parent=[&]() {
        int best=rng.index(n);
        for (int k=1;k<nt;++k) { int j=rng.index(n); if (values[j]<values[best]) best=j; }
        return pop[best];
    };
    for (int k=0;k<generations;++k) {
        int idx=static_cast<int>(std::min_element(values.begin(),values.end())-values.begin());
        Mat next(n,Vec(d)); Vec nv(n); next[0]=pop[idx]; nv[0]=values[idx];
        int count=1;
        while (count<n) {
            Vec a=parent(), b=parent();
            if (rng.u()<pc) for (std::size_t j=0;j<d;++j) if (rng.u()<0.5) std::swap(a[j],b[j]);
            for (Vec child : {a,b}) {
                if (count==n) break;
                if (rng.u()<pm) for (std::size_t j=0;j<d;++j) child[j]+=sigma[j]*rng.z();
                for (std::size_t j=0;j<d;++j) child[j]=clip(child[j],lo[j],hi[j]);
                next[count]=child; nv[count]=objective_value(f(child)); ++count;
            }
        }
        pop=next; values=nv;
    }
    int idx=static_cast<int>(std::min_element(values.begin(),values.end())-values.begin());
    Result r; r.point=pop[idx]; r.value=values[idx]; r.iterations=generations;
    r.evaluations=n+static_cast<long long>(generations)*(n-1); r.status="Search budget completed.";
    return r;
}

int main() {
  try {
    const unsigned int seed = 42;
    Random rng(seed);
    // Alternative: comment out the previous line and seed from system entropy.
    // Random rng(std::random_device{}());
    Function f=[](const Vec& p) { return std::sin(3*p[0])*std::cos(3*p[1])+(p[0]-1)*(p[0]-1)+(p[1]+2)*(p[1]+2); };
    Vec lower={-3,-3},upper={3,3},mutation_sigma={.24,.24};
    int population_size=100,max_generations=200,tournament_size=3;
    double crossover_probability=.9,mutation_probability=.3;
    print_result(genetic_search(f,lower,upper,rng,population_size,max_generations,crossover_probability,mutation_probability,mutation_sigma,tournament_size));
    std::cout << "Repeat with other seeds; global optimality is not guaranteed.\n";
    return 0;
  } catch (const std::exception& e) {
    std::cerr << "Error: " << e.what() << '\n'; return 1;
  }
}