nsga2.hpp

Go to the documentation of this file.

 
#pragma once
 
#include <algorithm>
#include <limits>
#include <numeric>
#include <stdexcept>
#include <vector>
 
#include "pulsatrix/individual.hpp"
 
namespace pulsatrix {
 
using Objectives = std::vector<double>;
 
inline bool Dominates(const Objectives& a, const Objectives& b) {
    if (a.size() != b.size()) {
        throw std::invalid_argument("Dominates: objective vectors must have equal size");
    }
    bool strictly_better_in_one = false;
    for (size_t i = 0; i < a.size(); ++i) {
        if (a[i] < b[i]) {
            return false;
        }
        if (a[i] > b[i]) {
            strictly_better_in_one = true;
        }
    }
    return strictly_better_in_one;
}
 
inline std::vector<std::vector<size_t>> FastNonDominatedSort(
    const std::vector<Objectives>& objectives) {
    const size_t n = objectives.size();
    std::vector<size_t> domination_count(n, 0);
    std::vector<std::vector<size_t>> dominates_set(n);
    std::vector<std::vector<size_t>> fronts;
    std::vector<size_t> current_front;
 
    for (size_t p = 0; p < n; ++p) {
        for (size_t q = 0; q < n; ++q) {
            if (p == q) {
                continue;
            }
            if (Dominates(objectives[p], objectives[q])) {
                dominates_set[p].push_back(q);
            } else if (Dominates(objectives[q], objectives[p])) {
                domination_count[p]++;
            }
        }
        if (domination_count[p] == 0) {
            current_front.push_back(p);
        }
    }
 
    while (!current_front.empty()) {
        fronts.push_back(current_front);
        std::vector<size_t> next_front;
        for (size_t p : current_front) {
            for (size_t q : dominates_set[p]) {
                if (--domination_count[q] == 0) {
                    next_front.push_back(q);
                }
            }
        }
        current_front = std::move(next_front);
    }
    return fronts;
}
 
inline std::vector<double> CrowdingDistance(const std::vector<Objectives>& front_objectives) {
    const size_t n = front_objectives.size();
    std::vector<double> distances(n, 0.0);
    if (n == 0) {
        return distances;
    }
    const size_t num_objectives = front_objectives[0].size();
    for (const auto& obj : front_objectives) {
        if (obj.size() != num_objectives) {
            throw std::invalid_argument(
                "CrowdingDistance: every objective vector in a front must have equal size");
        }
    }
    if (n <= 2) {
        std::fill(distances.begin(), distances.end(), std::numeric_limits<double>::infinity());
        return distances;
    }
 
    for (size_t m = 0; m < num_objectives; ++m) {
        std::vector<size_t> order(n);
        std::iota(order.begin(), order.end(), size_t{0});
        std::stable_sort(order.begin(), order.end(), [&](size_t a, size_t b) {
            return front_objectives[a][m] < front_objectives[b][m];
        });
 
        const double min_val = front_objectives[order.front()][m];
        const double max_val = front_objectives[order.back()][m];
        const double range = max_val - min_val;
 
        distances[order.front()] = std::numeric_limits<double>::infinity();
        distances[order.back()] = std::numeric_limits<double>::infinity();
 
        if (range <= 0.0) {
            continue;  // objective m is uniform across this front -- no information to add
        }
        for (size_t i = 1; i + 1 < n; ++i) {
            if (distances[order[i]] == std::numeric_limits<double>::infinity()) {
                continue;  // already a boundary individual under a different objective
            }
            distances[order[i]] +=
                (front_objectives[order[i + 1]][m] - front_objectives[order[i - 1]][m]) / range;
        }
    }
    return distances;
}
 
template <typename Genotype>
std::vector<Individual<Genotype, Objectives>> NSGA2Replacement(
    const std::vector<Individual<Genotype, Objectives>>& population,
    std::vector<Individual<Genotype, Objectives>> offspring, size_t mu) {
    if (population.size() + offspring.size() < mu) {
        throw std::invalid_argument(
            "NSGA2Replacement: population.size() + offspring.size() must be >= mu");
    }
 
    std::vector<Individual<Genotype, Objectives>> combined = population;
    combined.insert(combined.end(), std::make_move_iterator(offspring.begin()),
                     std::make_move_iterator(offspring.end()));
 
    std::vector<Objectives> objectives;
    objectives.reserve(combined.size());
    for (const auto& ind : combined) {
        objectives.push_back(ind.fitness);
    }
    auto fronts = FastNonDominatedSort(objectives);
 
    std::vector<Individual<Genotype, Objectives>> next_generation;
    next_generation.reserve(mu);
    for (const auto& front : fronts) {
        if (next_generation.size() + front.size() <= mu) {
            for (size_t idx : front) {
                next_generation.push_back(combined[idx]);
            }
            if (next_generation.size() == mu) {
                break;
            }
        } else {
            std::vector<Objectives> front_objectives;
            front_objectives.reserve(front.size());
            for (size_t idx : front) {
                front_objectives.push_back(combined[idx].fitness);
            }
            auto distances = CrowdingDistance(front_objectives);
 
            std::vector<size_t> order(front.size());
            std::iota(order.begin(), order.end(), size_t{0});
            std::stable_sort(order.begin(), order.end(),
                              [&](size_t a, size_t b) { return distances[a] > distances[b]; });
 
            const size_t remaining = mu - next_generation.size();
            for (size_t i = 0; i < remaining; ++i) {
                next_generation.push_back(combined[front[order[i]]]);
            }
            break;
        }
    }
    return next_generation;
}
 
}  // namespace pulsatrix