neat_evolution.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <cmath>
#include <limits>
#include <numeric>
#include <random>
#include <stdexcept>
#include <vector>
 
#include "pulsatrix/neat_genome.hpp"
#include "pulsatrix/neat_speciation.hpp"
 
namespace pulsatrix {
 
inline std::vector<int> AllocateOffspringCounts(const std::vector<double>& species_adjusted_fitness_sums,
                                                 int population_size) {
    if (species_adjusted_fitness_sums.empty()) {
        throw std::invalid_argument("AllocateOffspringCounts: species_adjusted_fitness_sums must not be empty");
    }
    if (population_size <= 0) {
        throw std::invalid_argument("AllocateOffspringCounts: population_size must be positive");
    }
 
    size_t n = species_adjusted_fitness_sums.size();
    double total = std::accumulate(species_adjusted_fitness_sums.begin(), species_adjusted_fitness_sums.end(), 0.0);
 
    std::vector<int> counts(n, 0);
    if (total <= 0.0) {
        int base = population_size / static_cast<int>(n);
        int remainder = population_size - base * static_cast<int>(n);
        for (size_t i = 0; i < n; ++i) {
            counts[i] = base + (static_cast<int>(i) < remainder ? 1 : 0);
        }
        return counts;
    }
 
    std::vector<double> exact(n);
    int allocated = 0;
    for (size_t i = 0; i < n; ++i) {
        exact[i] = species_adjusted_fitness_sums[i] / total * static_cast<double>(population_size);
        counts[i] = static_cast<int>(std::floor(exact[i]));
        allocated += counts[i];
    }
 
    int remaining = population_size - allocated;
    std::vector<size_t> order(n);
    std::iota(order.begin(), order.end(), 0);
    std::stable_sort(order.begin(), order.end(), [&](size_t a, size_t b) {
        return (exact[a] - std::floor(exact[a])) > (exact[b] - std::floor(exact[b]));
    });
    for (int i = 0; i < remaining; ++i) {
        counts[order[i]] += 1;
    }
    return counts;
}
 
template <typename RNG>
NEATGenome ReproduceOffspring(const NEATGenome& parent, InnovationTracker& tracker, RNG& rng,
                               double weight_mutation_sigma, double weight_mutation_probability,
                               double add_connection_probability, double add_node_probability) {
    NEATGenome child = parent;
    child.MutateWeights(weight_mutation_sigma, weight_mutation_probability, rng);
    std::bernoulli_distribution add_connection_roll(add_connection_probability);
    if (add_connection_roll(rng)) {
        child.AddConnection(tracker, rng);
    }
    std::bernoulli_distribution add_node_roll(add_node_probability);
    if (add_node_roll(rng)) {
        child.AddNode(tracker, rng);
    }
    return child;
}
 
struct NEATEvolutionResult {
    NEATGenome best_genome;
    double best_fitness;
    int generations_run;
};
 
template <typename FitnessFn, typename RNG>
NEATEvolutionResult RunNEATEvolution(std::vector<NEATGenome> population, FitnessFn fitness_fn, int num_generations,
                                      double compatibility_threshold, double c1, double c2, double c3,
                                      double weight_mutation_sigma, double weight_mutation_probability,
                                      double add_connection_probability, double add_node_probability,
                                      InnovationTracker& tracker, RNG& rng) {
    if (population.empty()) {
        throw std::invalid_argument("RunNEATEvolution: population must not be empty");
    }
    if (num_generations <= 0) {
        throw std::invalid_argument("RunNEATEvolution: num_generations must be positive");
    }
 
    NEATGenome best_genome = population[0];
    double best_fitness = -std::numeric_limits<double>::infinity();
 
    for (int gen = 0; gen < num_generations; ++gen) {
        std::vector<double> fitnesses(population.size());
        for (size_t i = 0; i < population.size(); ++i) {
            fitnesses[i] = fitness_fn(population[i]);
            if (fitnesses[i] > best_fitness) {
                best_fitness = fitnesses[i];
                best_genome = population[i];
            }
        }
 
        if (gen == num_generations - 1) {
            break;
        }
 
        SpeciesAssignment assignment = SpeciatePopulation(population, compatibility_threshold, c1, c2, c3);
        std::vector<double> adjusted = ComputeAdjustedFitness(fitnesses, assignment.species);
 
        std::vector<double> species_sums(assignment.species.size(), 0.0);
        for (size_t s = 0; s < assignment.species.size(); ++s) {
            for (size_t idx : assignment.species[s]) {
                species_sums[s] += adjusted[idx];
            }
        }
        std::vector<int> counts = AllocateOffspringCounts(species_sums, static_cast<int>(population.size()));
 
        std::vector<NEATGenome> next_population;
        next_population.reserve(population.size());
        for (size_t s = 0; s < assignment.species.size(); ++s) {
            int count = counts[s];
            if (count <= 0) {
                continue;
            }
            std::vector<size_t> members = assignment.species[s];
            std::sort(members.begin(), members.end(),
                      [&](size_t a, size_t b) { return fitnesses[a] > fitnesses[b]; });
 
            next_population.push_back(population[members.front()]);
            std::uniform_int_distribution<size_t> pick_parent(0, members.size() - 1);
            for (int offspring_index = 1; offspring_index < count; ++offspring_index) {
                const NEATGenome& parent = population[members[pick_parent(rng)]];
                next_population.push_back(ReproduceOffspring(parent, tracker, rng, weight_mutation_sigma,
                                                               weight_mutation_probability,
                                                               add_connection_probability, add_node_probability));
            }
        }
 
        population = std::move(next_population);
    }
 
    return NEATEvolutionResult{best_genome, best_fitness, num_generations};
}
 
}  // namespace pulsatrix