From 98ea1e2c0fc3eae20909fde19d758b392463d0b5 Mon Sep 17 00:00:00 2001 From: IWNMWE Date: Wed, 24 Jul 2024 16:54:02 +0530 Subject: [PATCH 01/11] added nsga3 declaration --- include/ensmallen_bits/nsga3/nsga3.hpp | 474 +++++++++++++++ include/ensmallen_bits/nsga3/nsga3_impl.hpp | 602 ++++++++++++++++++++ 2 files changed, 1076 insertions(+) create mode 100644 include/ensmallen_bits/nsga3/nsga3.hpp create mode 100644 include/ensmallen_bits/nsga3/nsga3_impl.hpp diff --git a/include/ensmallen_bits/nsga3/nsga3.hpp b/include/ensmallen_bits/nsga3/nsga3.hpp new file mode 100644 index 000000000..0d27dadc5 --- /dev/null +++ b/include/ensmallen_bits/nsga3/nsga3.hpp @@ -0,0 +1,474 @@ +/** + * @file nsga3.hpp + * @author Satyam Shukla + * + * NSGA-III is a multi-objective optimization algorithm, widely used in + * many real-world applications. NSGA-III generates offsprings using + * crossover and mutation and then selects the next generation according + * to non-dominated-sorting and crowding distance comparison.The maintenance + * of diversity among population members in NSGA-III is aided by supplying and + * adaptively updating a number of well-spread reference points. + * + * ensmallen is free software; you may redistribute it and/or modify it under + * the terms of the 3-clause BSD license. You should have received a copy of + * the 3-clause BSD license along with ensmallen. If not, see + * http://www.opensource.org/licenses/BSD-3-Clause for more information. + */ + +#ifndef ENSMALLEN_NSGA3_NSGA3_HPP +#define ENSMALLEN_NSGA3_NSGA3_HPP + +namespace ens { + +/** + * NSGA-III (Non-dominated Sorting Genetic Algorithm - III) is a multi-objective + * optimization algorithm. This class implements the NSGA-III algorithm. + * + * The algorithm works by generating a candidate population from a fixed + * starting point. At each stage of optimization, a new population of children + * is generated. This new population along with its predecessor is sorted using + * non-domination as the metric. Following this, the population is further + * segregated in fronts. A new population is generated from these fronts having + * size equal to that of the starting population. + * + * During evolution, two parents are randomly chosen using binary tournament + * selection. A pair of children are generated by crossing over these two + * candidates followed by mutation. + * + * The best front (Pareto optimal) is returned by the Optimize() method. + * + * For more information, see the following: + * + * @code + * @article{deb2013evolutionary, + * title={An evolutionary many-objective optimization algorithm using reference-point-based nondominated sorting approach, part I: solving problems with box constraints}, + * author={Deb, Kalyanmoy and Jain, Himanshu}, + * journal={IEEE transactions on evolutionary computation}, + * volume={18}, + * number={4}, + * pages={577--601}, + * year={2013}, + * publisher={IEEE} + * } + * @endcode + * + * NSGA-III can optimize arbitrary multi-objective functions. For more details, + * see the documentation on function types included with this distribution or + * on the ensmallen website. + */ +template +class NSGA3 +{ + public: + /** + * Constructor for the NSGA3 optimizer. + * + * The default values provided over here are not necessarily suitable for a + * given function. Therefore it is highly recommended to adjust the + * parameters according to the problem. + * + * @param referencePoints The reference points to be used. + * @param populationSize The number of candidates in the population. + * This should be atleast 4 in size and a multiple of 4. + * @param maxGenerations The maximum number of generations allowed for NSGA-II. + * @param crossoverProb The probability that a crossover will occur. + * @param mutationProb The probability that a mutation will occur. + * @param mutationStrength The strength of the mutation. + * @param epsilon The minimum difference required to distinguish between + * candidate solutions. + * @param lowerBound Lower bound of the coordinates of the initial population. + * @param upperBound Upper bound of the coordinates of the initial population. + */ + NSGA3(const MatType& referencePoints, + const size_t populationSize = 100, + const size_t maxGenerations = 2000, + const double crossoverProb = 0.6, + const double mutationProb = 0.3, + const double mutationStrength = 1e-3, + const double epsilon = 1e-6, + const arma::vec& lowerBound = arma::zeros(1, 1), + const arma::vec& upperBound = arma::ones(1, 1)); + + /** + * Constructor for the NSGA3 optimizer. This constructor provides an overload + * to use `lowerBound` and `upperBound` of type double. + * + * The default values provided over here are not necessarily suitable for a + * given function. Therefore it is highly recommended to adjust the + * parameters according to the problem. + * + * @param referencePoints The reference points to be used. + * @param populationSize The number of candidates in the population. + * This should be atleast 4 in size and a multiple of 4. + * @param maxGenerations The maximum number of generations allowed for NSGA-II. + * @param crossoverProb The probability that a crossover will occur. + * @param mutationProb The probability that a mutation will occur. + * @param mutationStrength The strength of the mutation. + * @param epsilon The minimum difference required to distinguish between + * candidate solutions. + * @param lowerBound Lower bound of the coordinates of the initial population. + * @param upperBound Upper bound of the coordinates of the initial population. + */ + NSGA3(const MatType& referencePoints, + const size_t populationSize = 100, + const size_t maxGenerations = 2000, + const double crossoverProb = 0.6, + const double mutationProb = 0.3, + const double mutationStrength = 1e-3, + const double epsilon = 1e-6, + const double lowerBound = 0, + const double upperBound = 1); + + /** + * Optimize a set of objectives. The initial population is generated using the + * starting point. The output is the best generated front. + * + * @tparam ArbitraryFunctionType std::tuple of multiple objectives. + * @tparam CallbackTypes Types of callback functions. + * @param objectives Vector of objective functions to optimize for. + * @param iterate Starting point. + * @param callbacks Callback functions. + * @return MatType::elem_type The minimum of the accumulated sum over the + * objective values in the best front. + */ + template + typename MatType::elem_type Optimize( + std::tuple& objectives, + MatType& iterate, + CallbackTypes&&... callbacks); + + //! Get the population size. + size_t PopulationSize() const { return populationSize; } + //! Modify the population size. + size_t& PopulationSize() { return populationSize; } + + //! Get the maximum number of generations. + size_t MaxGenerations() const { return maxGenerations; } + //! Modify the maximum number of generations. + size_t& MaxGenerations() { return maxGenerations; } + + //! Get the crossover rate. + double CrossoverRate() const { return crossoverProb; } + //! Modify the crossover rate. + double& CrossoverRate() { return crossoverProb; } + + //! Get the mutation probability. + double MutationProbability() const { return mutationProb; } + //! Modify the mutation probability. + double& MutationProbability() { return mutationProb; } + + //! Get the mutation strength. + double MutationStrength() const { return mutationStrength; } + //! Modify the mutation strength. + double& MutationStrength() { return mutationStrength; } + + //! Get the tolerance. + double Epsilon() const { return epsilon; } + //! Modify the tolerance. + double& Epsilon() { return epsilon; } + + //! Get the reference points. + const MatType& ReferencePoints() const { return referencePoints; } + //! Modify the reference points. + MatType& ReferencePoints() { return referencePoints; } + + //! Retrieve value of lowerBound. + const arma::vec& LowerBound() const { return lowerBound; } + //! Modify value of lowerBound. + arma::vec& LowerBound() { return lowerBound; } + + //! Retrieve value of upperBound. + const arma::vec& UpperBound() const { return upperBound; } + //! Modify value of upperBound. + arma::vec& UpperBound() { return upperBound; } + + //! Retrieve the Pareto optimal points in variable space. This returns an empty cube + //! until `Optimize()` has been called. + const arma::cube& ParetoSet() const { return paretoSet; } + + //! Retrieve the best front (the Pareto frontier). This returns an empty cube until + //! `Optimize()` has been called. + const arma::cube& ParetoFront() const { return paretoFront; } + + /** + * Retrieve the best front (the Pareto frontier). This returns an empty + * vector until `Optimize()` has been called. Note that this function is + * deprecated and will be removed in ensmallen 3.x! Use `ParetoFront()` + * instead. + */ + [[deprecated("use ParetoFront() instead")]] const std::vector& Front() + { + if (rcFront.size() == 0) + { + // Match the old return format. + for (size_t i = 0; i < paretoFront.n_slices; ++i) + { + rcFront.push_back(arma::mat(paretoFront.slice(i))); + } + } + + return rcFront; + } + + private: + /** + * Evaluate objectives for the elite population. + * + * @tparam ArbitraryFunctionType std::tuple of multiple function types. + * @tparam MatType Type of matrix to optimize. + * @param population The elite population. + * @param objectives The set of objectives. + * @param calculatedObjectives Vector to store calculated objectives. + */ + template + typename std::enable_if::type + EvaluateObjectives(std::vector&, + std::tuple&, + std::vector >&); + + template + typename std::enable_if::type + EvaluateObjectives(std::vector& population, + std::tuple& objectives, + std::vector >& + calculatedObjectives); + + /** + * Reproduce candidates from the elite population to generate a new + * population. + * + * @tparam MatType Type of matrix to optimize. + * @param population The elite population. + * @param objectives The set of objectives. + * @param lowerBound Lower bound of the coordinates of the initial population. + * @param upperBound Upper bound of the coordinates of the initial population. + */ + template + void BinaryTournamentSelection(std::vector& population, + const MatType& lowerBound, + const MatType& upperBound); + + /** + * Crossover two parents to create a pair of new children. + * + * @tparam MatType Type of matrix to optimize. + * @param childA A newly generated candidate. + * @param childB Another newly generated candidate. + * @param parentA First parent from elite population. + * @param parentB Second parent from elite population. + */ + template + void Crossover(MatType& childA, + MatType& childB, + const MatType& parentA, + const MatType& parentB); + + /** + * Mutate the coordinates for a candidate. + * + * @tparam MatType Type of matrix to optimize. + * @param child The candidate whose coordinates are being modified. + * @param objectives The set of objectives. + * @param lowerBound Lower bound of the coordinates of the initial population. + * @param upperBound Upper bound of the coordinates of the initial population. + */ + template + void Mutate(MatType& child, + const MatType& lowerBound, + const MatType& upperBound); + + /** + * Sort the candidate population using their domination count and the set of + * dominated nodes. + * + * @tparam MatType Type of matrix to optimize. + * @param fronts The population is sorted into these Pareto fronts. The first + * front is the best, the second worse and so on. + * @param ranks The assigned ranks, used for crowding distance based sorting. + * @param calculatedObjectives The previously calculated objectives. + */ + template + void FastNonDominatedSort( + std::vector >& fronts, + std::vector& ranks, + std::vector >& calculatedObjectives); + + /** + * Normalizes the front given the extreme points in the current front. + * + * @tparam The type of population datapoints. + * @param calculatedObjectives The current population evaluated objectives. + * @param normalization The normalizing vector. + * @param front The previously generated Pareto front. + * @param extreme The indexes of the extreme points in the front. + */ + void NormalizeFront( + std::vector >& calculatedObjectives, + arma::Col& normalization, + const std::vector& front, + const arma::Row& extreme); + + /** + * Finding the indexes of the extreme points in the front. + * + * @param indexes vector containing the slected indexes. + * @param calculatedObjectives The current population objectives. + * @param front The front of the current generation. + */ + void FindExtremePoints(arma::Row& indexes, + std::vector >& calculatedObjectives, + const std::vector& front); + + /** + * Finding the distance of each point in the front from the line formed + * by pointA and pointB. + * + * @param distance The vector containing the distances of the points in the fron from the line. + * @param calculatedObjectives Reference to the current population evaluated Objectives. + * @param front The front of the current generation(indices of population). + * @param pointA The first point on the line. + * @param pointB The second point on the line. + */ + void PointToLineDistance(arma::Row& distances, + std::vector >& calculatedObjectives, + const std::vector& front, + const arma::Col& pointA, + const arma::Col& pointB); + + /** + * Finding the point in the reference front associated with the members in + * the given front and also stores the distance between the two points. + * + * @param index Vector containing the index of the point in the refrence directons associated. + * @param dists Vector of distances from the corresponding point in the front to the associated reference direction. + * @param population The points of the currently generated population. + * @param front The index of points belonging to the given front. + */ + void Associate(arma::Row& refindex, + arma::Row& dists, + const std::vector& population, + const std::vector& front); + + /** + * Find the niche count for each reference direction. + * + * @param count The no. of points associated with each niche direction. + * @param refIndex Vector containing the index of the point in the reference + * directons associated. + */ + void NicheCount(arma::Row& count, + const arma::Row& refIndex); + + /** + * The niche preserving operation to select the final points form the given front + * aranges the front in descending order of priority for the top K points in the front. + * + * @param K The no. of remaining points to select from the given front for the next population. + * @param refIndex The index of the rerference points associated with the in the given front. + * @param dists The distances of th points in the front form their associated reference points line. + * @param front The index of teh points within the population which are a part of the given front. + */ + void Niching(size_t K. + const arma::Row& refindex, + const arma::Row& dists, + const std::vector& front); + + /** + * Operator to check if one candidate Pareto-dominates the other. + * + * A candidate is said to dominate the other if it is at least as good as the + * other candidate for all the objectives and there exists at least one + * objective for which it is strictly better than the other candidate. + * + * @tparam MatType Type of matrix to optimize. + * @param calculatedObjectives The previously calculated objectives. + * @param candidateP The candidate being compared from the elite population. + * @param candidateQ The candidate being compared against. + * @return true if candidateP Pareto dominates candidateQ, otherwise, false. + */ + template + bool Dominates( + std::vector >& calculatedObjectives, + size_t candidateP, + size_t candidateQ); + + /** + * The operator used in the crowding distance based sorting. + * + * If a candidates has a lower rank then it is preferred. + * Otherwise, if the ranks are equal then the candidate with the larger + * crowding distance is preferred. + * + * @param idxP The index of the first cadidate from the elite population being + * sorted. + * @param idxQ The index of the second cadidate from the elite population + * being sorted. + * @param ranks The previously calculated ranks. + * @param crowdingDistance The crowding distance for each individual in + * the population. + * @return true if the first candidate is preferred, otherwise, false. + */ + template + bool CrowdingOperator(size_t idxP, + size_t idxQ, + const std::vector& ranks, + const std::vector& crowdingDistance); + + //! The number of objectives being optimised for. + size_t numObjectives; + + //! The numbeer of variables used per objectives. + size_t numVariables; + + //! The number of candidates in the population. + size_t populationSize; + + //! Maximum number of generations before termination criteria is met. + size_t maxGenerations; + + //! Probability that crossover will occur. + double crossoverProb; + + //! Probability that mutation will occur. + double mutationProb; + + //! Strength of the mutation. + double mutationStrength; + + //! The tolerance for termination. + double epsilon; + + //! Lower bound of the initial swarm. + arma::vec lowerBound; + + //! Upper bound of the initial swarm. + arma::vec upperBound; + + //! The set of all the Pareto optimal points. + //! Stored after Optimize() is called. + arma::cube paretoSet; + + //! The set of all the Pareto optimal objective vectors. + //! Stored after Optimize() is called. + arma::cube paretoFront; + + //! The reference points. + MatType referencePoints; + + //! A different representation of the Pareto front, for reverse compatibility + //! purposes. This can be removed when ensmallen 3.x is released! (Along + //! with `Front()`.) This is only populated when `Front()` is called. + std::vector rcFront; +}; + +} // namespace ens + +// Include implementation. +#include "nsga3_impl.hpp" + +#endif diff --git a/include/ensmallen_bits/nsga3/nsga3_impl.hpp b/include/ensmallen_bits/nsga3/nsga3_impl.hpp new file mode 100644 index 000000000..7cb04fdfc --- /dev/null +++ b/include/ensmallen_bits/nsga3/nsga3_impl.hpp @@ -0,0 +1,602 @@ +/** + * @file nsga3_impl.hpp + * @author Satyam Shukla + * + * Implementation of the NSGA3 algorithm. Used for multi-objective + * optimization problems on arbitrary functions. + * + * ensmallen is free software; you may redistribute it and/or modify it under + * the terms of the 3-clause BSD license. You should have received a copy of + * the 3-clause BSD license along with ensmallen. If not, see + * http://www.opensource.org/licenses/BSD-3-Clause for more Information. + */ + +#ifndef ENSMALLEN_NSGA3_NSGA3_IMPL_HPP +#define ENSMALLEN_NSGA3_NSGA3_IMPL_HPP + +#include "nsga3.hpp" +#include + +namespace ens { + +template +inline NSGA3::NSGA3(const MatType& referencePoints, + const size_t populationSize = 100, + const size_t maxGenerations = 2000, + const double crossoverProb = 0.6, + const double mutationProb = 0.3, + const double mutationStrength = 1e-3, + const double epsilon = 1e-6, + const arma::vec& lowerBound = arma::zeros(1, 1), + const arma::vec& upperBound = arma::ones(1, 1)): + referencePoints(referencePoints), + numObjectives(0), + numVariables(0), + populationSize(populationSize), + maxGenerations(maxGenerations), + crossoverProb(crossoverProb), + distributionIndex(distributionIndex), + epsilon(epsilon), + eta(eta), + lowerBound(lowerBound), + upperBound(upperBound) +{ /* Nothing to do here. */ } + +template +inline NSGA3::NSGA3(const MatType& referencePoints, + const size_t populationSize = 100, + const size_t maxGenerations = 2000, + const double crossoverProb = 0.6, + const double mutationProb = 0.3, + const double mutationStrength = 1e-3, + const double epsilon = 1e-6, + const double lowerBound = 0, + const double upperBound = 1): + referencePoints(referencePoints), + numObjectives(0), + numVariables(0), + populationSize(populationSize), + maxGenerations(maxGenerations), + crossoverProb(crossoverProb), + distributionIndex(distributionIndex), + epsilon(epsilon), + eta(eta), + lowerBound(lowerBound * arma::ones(1, 1)), + upperBound(upperBound * arma::ones(1, 1)) +{ /* Nothing to do here. */ } + +//! Optimize the function. +template +template +typename MatType::elem_type NSGA3::Optimize( + std::tuple& objectives, + MatType& iterateIn, + CallbackTypes&&... callbacks) +{ + // Make sure for evolution to work at least four candidates are present. + if (populationSize < 4 && populationSize % 4 != 0) + { + throw std::logic_error("NSGA3::Optimize(): population size should be at" + " least 4, and, a multiple of 4!"); + } + + // Convenience typedefs. + typedef typename MatType::elem_type ElemType; + typedef typename MatTypeTraits::BaseMatType BaseMatType; + + BaseMatType& iterate = (BaseMatType&) iterateIn; + + // Make sure that we have the methods that we need. Long name... + traits::CheckArbitraryFunctionTypeAPI(); + RequireDenseFloatingPointType(); + + // Check if lower bound is a vector of a single dimension. + if (lowerBound.n_rows == 1) + lowerBound = lowerBound(0, 0) * arma::ones(iterate.n_rows, iterate.n_cols); + + // Check if upper bound is a vector of a single dimension. + if (upperBound.n_rows == 1) + upperBound = upperBound(0, 0) * arma::ones(iterate.n_rows, iterate.n_cols); + + // Check the dimensions of lowerBound and upperBound. + assert(lowerBound.n_rows == iterate.n_rows && "The dimensions of " + "lowerBound are not the same as the dimensions of iterate."); + assert(upperBound.n_rows == iterate.n_rows && "The dimensions of " + "upperBound are not the same as the dimensions of iterate."); + + numObjectives = sizeof...(ArbitraryFunctionType); + numVariables = iterate.n_rows; + + // Cache calculated objectives. + std::vector > calculatedObjectives(populationSize); + + // Population size reserved to 2 * populationSize + 1 to accommodate + // for the size of intermediate candidate population. + std::vector population; + population.reserve(2 * populationSize + 1); + + // Pareto fronts, initialized during non-dominated sorting. + // Stores indices of population belonging to a certain front. + std::vector > fronts; + // Initialised in SurvivalScoreAssignment. + std::vector survivalScore; + // Initialised during non-dominated sorting. + std::vector ranks; + + //! Useful temporaries for float-like comparisons. + const BaseMatType castedLowerBound = arma::conv_to::from(lowerBound); + const BaseMatType castedUpperBound = arma::conv_to::from(upperBound); + + // Controls early termination of the optimization process. + bool terminate = false; + + // Generate the population based on a uniform distribution around the given + // starting point. + for (size_t i = 0; i < populationSize; i++) + { + population.push_back(arma::randu(iterate.n_rows, + iterate.n_cols) - 0.5 + iterate); + + // Constrain all genes to be within bounds. + population[i] = arma::min(arma::max(population[i], castedLowerBound), castedUpperBound); + } + + Info << "NSGA3 initialized successfully. Optimization started." << std::endl; + + // Iterate until maximum number of generations is obtained. + Callback::BeginOptimization(*this, objectives, iterate, callbacks...); + + for (size_t generation = 1; generation <= maxGenerations && !terminate; generation++) + { + // Create new population of candidate from the present elite population. + // Have P_t, generate G_t using P_t. + BinaryTournamentSelection(population, castedLowerBound, castedUpperBound); + + // Evaluate the objectives for the new population. + calculatedObjectives.resize(population.size()); + std::fill(calculatedObjectives.begin(), calculatedObjectives.end(), + arma::Col(numObjectives, arma::fill::zeros)); + EvaluateObjectives(population, objectives, calculatedObjectives); + + // Perform fast non dominated sort on P_t ∪ G_t. + ranks.resize(population.size()); + FastNonDominatedSort(fronts, ranks, calculatedObjectives); + + std::vector selectedPoints(fronts[0].begin(), fronts[0].end()); + while (selectedPoints.size() < populationSize) + { + selectedPoints.insert(selectedPoints.end(), fronts[index].begin(), fronts[index].end()); + } + arma::Col idealPoint(calculatedObjectives[selectedPoints[0]]); + for (size_t index = 1; index < selectedPoints.size(); index++) + { + idealPoint = arma::min(idealPoint, calculatedObjectives[selectedPoints[index]]); + } + + arma::Col normalize(numObjectives, arma::fill::zeros); + + + // Sort based on survival score. + std::sort(population.begin(), population.end(), + [this, ranks, survivalScore, population] + (BaseMatType candidateP, BaseMatType candidateQ) + { + size_t idxP{}, idxQ{}; + for (size_t i = 0; i < population.size(); i++) + { + if (arma::approx_equal(population[i], candidateP, "absdiff", epsilon)) + idxP = i; + + if (arma::approx_equal(population[i], candidateQ, "absdiff", epsilon)) + idxQ = i; + } + + return SurvivalScoreOperator(idxP, idxQ, ranks, survivalScore); + } + ); + + // Yield a new population P_{t+1} of size populationSize. + // Discards unfit population from the R_{t} to yield P_{t+1}. + population.resize(populationSize); + + terminate |= Callback::GenerationalStepTaken(*this, objectives, iterate, + calculatedObjectives, fronts, callbacks...); + } + EvaluateObjectives(population, objectives, calculatedObjectives); + // Set the candidates from the Pareto Set as the output. + paretoSet.set_size(population[0].n_rows, population[0].n_cols, fronts[0].size()); + // The Pareto Set is stored, can be obtained via ParetoSet() getter. + for (size_t solutionIdx = 0; solutionIdx < fronts[0].size(); ++solutionIdx) + { + paretoSet.slice(solutionIdx) = + arma::conv_to::from(population[fronts[0][solutionIdx]]); + } + + // Set the candidates from the Pareto Front as the output. + paretoFront.set_size(calculatedObjectives[0].n_rows, calculatedObjectives[0].n_cols, + fronts[0].size()); + // The Pareto Front is stored, can be obtained via ParetoFront() getter. + for (size_t solutionIdx = 0; solutionIdx < fronts[0].size(); ++solutionIdx) + { + paretoFront.slice(solutionIdx) = + arma::conv_to::from(calculatedObjectives[fronts[0][solutionIdx]]); + } + + // Clear rcFront, in case it is later requested by the user for reverse + // compatibility reasons. + rcFront.clear(); + + // Assign iterate to first element of the Pareto Set. + iterate = population[fronts[0][0]]; + + Callback::EndOptimization(*this, objectives, iterate, callbacks...); + + ElemType performance = std::numeric_limits::max(); + + for (const arma::Col& objective: calculatedObjectives) + if (arma::accu(objective) < performance) + performance = arma::accu(objective); + + return performance; +} + +//! No objectives to evaluate. +template +typename std::enable_if::type +AGEMOEA::EvaluateObjectives( + std::vector&, + std::tuple&, + std::vector >&) +{ + // Nothing to do here. +} + +//! Evaluate the objectives for the entire population. +template +typename std::enable_if::type +AGEMOEA::EvaluateObjectives( + std::vector& population, + std::tuple& objectives, + std::vector >& calculatedObjectives) +{ + for (size_t i = 0; i < populationSize; i++) + { + calculatedObjectives[i](I) = std::get(objectives).Evaluate(population[i]); + EvaluateObjectives(population, objectives, + calculatedObjectives); + } +} + +//! Reproduce and generate new candidates. +template +inline void AGEMOEA::BinaryTournamentSelection(std::vector& population, + const MatType& lowerBound, + const MatType& upperBound) +{ + std::vector children; + + while (children.size() < population.size()) + { + // Choose two random parents for reproduction from the elite population. + size_t indexA = arma::randi(arma::distr_param(0, populationSize - 1)); + size_t indexB = arma::randi(arma::distr_param(0, populationSize - 1)); + + // Make sure that the parents differ. + if (indexA == indexB) + { + if (indexB < populationSize - 1) + indexB++; + else + indexB--; + } + + // Initialize the children to the respective parents. + MatType childA = population[indexA], childB = population[indexB]; + + if(arma::randu() <= crossoverProb) + Crossover(childA, childB, population[indexA], population[indexB], + lowerBound, upperBound); + + Mutate(childA, 1.0 / static_cast(numVariables), + lowerBound, upperBound); + Mutate(childB, 1.0 / static_cast(numVariables), + lowerBound, upperBound); + + // Add the children to the candidate population. + children.push_back(childA); + children.push_back(childB); + } + + // Add the candidates to the elite population. + population.insert(std::end(population), std::begin(children), std::end(children)); +} + +//! Perform simulated binary crossover (SBX) of genes for the children. +template +inline void AGEMOEA::Crossover(MatType& childA, + MatType& childB, + const MatType& parentA, + const MatType& parentB, + const MatType& lowerBound, + const MatType& upperBound) +{ + //! Generates a child from two parent individuals + // according to the polynomial probability distribution. + arma::Cube parents(parentA.n_rows, parentA.n_cols, 2); + parents.slice(0) = parentA; + parents.slice(1) = parentB; + MatType current_min = arma::min(parents, 2); + MatType current_max = arma::max(parents, 2); + + if (arma::accu(parentA - parentB < 1e-14)) + { + childA = parentA; + childB = parentB; + return; + } + MatType current_diff = current_max - current_min; + current_diff.transform( [](typename MatType::elem_type val) + { return (val < 1e-10 ? 1e-10:val); } ); + + // Calculating beta used for the final crossover. + MatType beta1 = 1 + 2.0 * (current_min - lowerBound) / current_diff; + MatType beta2 = 1 + 2.0 * (upperBound - current_max) / current_diff; + MatType alpha1 = 2 - arma::pow(beta1, -(eta + 1)); + MatType alpha2 = 2 - arma::pow(beta2, -(eta + 1)); + + MatType us(arma::size(alpha1), arma::fill::randu); + arma::umat mask1 = us > (1.0 / alpha1); + MatType betaq1 = arma::pow(us % alpha1, 1. / (eta + 1)); + betaq1 = betaq1 % (mask1 != 1.0) + arma::pow((1.0 / (2.0 - us % alpha1)), 1.0 / (eta + 1)) % mask1; + arma::umat mask2 = us > (1.0 / alpha2); + MatType betaq2 = arma::pow(us % alpha2, 1 / (eta + 1)); + betaq2 = betaq2 % (mask1 != 1.0) + arma::pow((1.0 / (2.0 - us % alpha2)), 1.0 / (eta + 1)) % mask2; + + // Variables after the cross over for all of them. + MatType c1 = 0.5 * ((current_min + current_max) - betaq1 % current_diff); + MatType c2 = 0.5 * ((current_min + current_max) + betaq2 % current_diff); + c1 = arma::min(arma::max(c1, lowerBound), upperBound); + c2 = arma::min(arma::max(c2, lowerBound), upperBound); + + // Decision for the crossover between the two parents for each variable. + us.randu(); + childA = parentA % (us <= 0.5); + childB = parentB % (us <= 0.5); + us.randu(); + childA = childA + c1 % ((us <= 0.5) % (childA == 0)); + childA = childA + c2 % ((us > 0.5) % (childA == 0)); + childB = childB + c2 % ((us <= 0.5) % (childB == 0)); + childB = childB + c1 % ((us > 0.5) % (childB == 0)); +} + +//! Perform Polynomial mutation of the candidate. +template +inline void AGEMOEA::Mutate(MatType& candidate, + double mutationRate, + const MatType& lowerBound, + const MatType& upperBound) +{ + const size_t numVariables = candidate.n_rows; + for (size_t geneIdx = 0; geneIdx < numVariables; ++geneIdx) + { + // Should this gene be mutated? + if (arma::randu() > mutationRate) + continue; + + const double geneRange = upperBound(geneIdx) - lowerBound(geneIdx); + // Normalised distance from the bounds. + const double lowerDelta = (candidate(geneIdx) - lowerBound(geneIdx)) / geneRange; + const double upperDelta = (upperBound(geneIdx) - candidate(geneIdx)) / geneRange; + const double mutationPower = 1. / (distributionIndex + 1.0); + const double rand = arma::randu(); + double value, perturbationFactor; + if (rand < 0.5) + { + value = 2.0 * rand + (1.0 - 2.0 * rand) * + std::pow(upperDelta, distributionIndex + 1.0); + perturbationFactor = std::pow(value, mutationPower) - 1.0; + } + else + { + value = 2.0 * (1.0 - rand) + 2.0 *(rand - 0.5) * + std::pow(lowerDelta, distributionIndex + 1.0); + perturbationFactor = 1.0 - std::pow(value, mutationPower); + } + + candidate(geneIdx) += perturbationFactor * geneRange; + } + //! Enforce bounds. + candidate = arma::min(arma::max(candidate, lowerBound), upperBound); +} + +template +inline void AGEMOEA::NormalizeFront( + std::vector >& calculatedObjectives, + arma::Col& normalization, + const std::vector& front, + const arma::Row& extreme) +{ + arma::Mat vectorizedObjectives(numObjectives, front.size()); + for (size_t i = 0; i < front.size(); i++) + { + vectorizedObjectives.col(i) = calculatedObjectives[front[i]]; + } + + if(front.size() < numObjectives) + { + normalization = arma::max(vectorizedObjectives, 1); + return; + } + arma::Col temp; + arma::uvec unique = arma::find_unique(extreme); + if (extreme.n_elem != unique.n_elem) + { + normalization = arma::max(vectorizedObjectives, 1); + return; + } + arma::Col one(front.size(), arma::fill::ones); + arma::Col hyperplane = arma::solve( + vectorizedObjectives.t(), one); + if (hyperplane.has_inf() || hyperplane.has_nan() || (arma::accu(hyperplane < 0.0) > 0)) + { + normalization = arma::max(vectorizedObjectives, 1); + } + else + { + normalization = 1. / hyperplane; + if (hyperplane.has_inf() || hyperplane.has_nan()) + { + normalization = arma::max(vectorizedObjectives, 1); + } + } + normalization = normalization + (normalization == 0); +} + +//! Find the index of the of the extreme points in the given front. +template +void NSGA::FindExtremePoints(arma::Row& indexes, + std::vector >& calculatedObjectives, + const std::vector& front) +{ + typedef typename MatType::elem_type ElemType; + + if(numObjectives >= front.size()) + { + indexes = arma::linspace>(0, front.size() - 1, front.size()); + return; + } + + arma::Mat W(numObjectives, numObjectives, arma::fill::eye); + W = W + 1e-6; + std::vector selected(front.size()); + arma::Col z(numObjectives, arma::fill::zeros); + arma::Row dists; + for (size_t i = 0; i < numObjectives; i++) + { + PointToLineDistance(dists, calculatedObjectives, front, z, W.col(i)); + for (size_t j = 0; j < front.size(); j++) + if (selected[j]){dists[j] = arma::datum::inf;} + indexes[i] = dists.index_min(); + selected[dists.index_min()] = true; + } +} + +//! Find the distance of a front from a line formed by two points. +template +void NSGA3::PointToLineDistance(arma::Row& distances, + std::vector >& calculatedObjectives, + const std::vector& front, + const arma::Col& pointA, + const arma::Col& pointB) +{ + typedef typename MatType::elem_type ElemType; + arma::Row distancesTemp(front.size()); + arma::Col ba = pointB - pointA; + arma::Col pa; + + for (size_t i = 0; i < front.size(); i++) + { + size_t ind = front[i]; + + pa = (calculatedObjectives[ind] - pointA); + double t = arma::dot(pa, ba) / arma::dot(ba, ba); + distancesTemp[i] = arma::accu(arma::pow((pa - t * ba), 2)); + } + distances = distancesTemp; +} + +//! Sort population into Pareto fronts. +template +inline void AGEMOEA::FastNonDominatedSort( + std::vector >& fronts, + std::vector& ranks, + std::vector >& calculatedObjectives) +{ + std::map dominationCount; + std::map > dominated; + + // Reset and initialize fronts. + fronts.clear(); + fronts.push_back(std::vector()); + + for (size_t p = 0; p < populationSize; p++) + { + dominated[p] = std::set(); + dominationCount[p] = 0; + + for (size_t q = 0; q < populationSize; q++) + { + if (Dominates(calculatedObjectives, p, q)) + dominated[p].insert(q); + else if (Dominates(calculatedObjectives, q, p)) + dominationCount[p] += 1; + } + + if (dominationCount[p] == 0) + { + ranks[p] = 0; + fronts[0].push_back(p); + } + } + + size_t i = 0; + + while (!fronts[i].empty()) + { + std::vector nextFront; + + for (size_t p: fronts[i]) + { + for (size_t q: dominated[p]) + { + dominationCount[q]--; + + if (dominationCount[q] == 0) + { + ranks[q] = i + 1; + nextFront.push_back(q); + } + } + } + + i++; + fronts.push_back(nextFront); + } + // Remove the empty final set. + fronts.pop_back(); +} + +//! Check if a candidate Pareto dominates another candidate. +template +inline bool AGEMOEA::Dominates( + std::vector >& calculatedObjectives, + size_t candidateP, + size_t candidateQ) +{ + bool allBetterOrEqual = true; + bool atleastOneBetter = false; + size_t n_objectives = calculatedObjectives[0].n_elem; + + for (size_t i = 0; i < n_objectives; i++) + { + // P is worse than Q for the i-th objective function. + if (calculatedObjectives[candidateP](i) > calculatedObjectives[candidateQ](i)) + allBetterOrEqual = false; + + // P is better than Q for the i-th objective function. + else if (calculatedObjectives[candidateP](i) < calculatedObjectives[candidateQ](i)) + atleastOneBetter = true; + } + + return allBetterOrEqual && atleastOneBetter; +} + +} // namespace ens + +#endif \ No newline at end of file From e5b4d3a1717267e4b904f9a460b58f3c04ddc0b0 Mon Sep 17 00:00:00 2001 From: IWNMWE Date: Fri, 2 Aug 2024 16:44:30 +0530 Subject: [PATCH 02/11] nsga3 definition --- include/ensmallen_bits/nsga3/nsga3.hpp | 97 +++---- include/ensmallen_bits/nsga3/nsga3_impl.hpp | 281 +++++++++++++------- 2 files changed, 227 insertions(+), 151 deletions(-) diff --git a/include/ensmallen_bits/nsga3/nsga3.hpp b/include/ensmallen_bits/nsga3/nsga3.hpp index 0d27dadc5..0831e9257 100644 --- a/include/ensmallen_bits/nsga3/nsga3.hpp +++ b/include/ensmallen_bits/nsga3/nsga3.hpp @@ -79,7 +79,7 @@ class NSGA3 * @param lowerBound Lower bound of the coordinates of the initial population. * @param upperBound Upper bound of the coordinates of the initial population. */ - NSGA3(const MatType& referencePoints, + NSGA3(const arma::Cube& referencePoints, const size_t populationSize = 100, const size_t maxGenerations = 2000, const double crossoverProb = 0.6, @@ -109,7 +109,7 @@ class NSGA3 * @param lowerBound Lower bound of the coordinates of the initial population. * @param upperBound Upper bound of the coordinates of the initial population. */ - NSGA3(const MatType& referencePoints, + NSGA3(const arma::Cube& referencePoints, const size_t populationSize = 100, const size_t maxGenerations = 2000, const double crossoverProb = 0.6, @@ -169,9 +169,9 @@ class NSGA3 double& Epsilon() { return epsilon; } //! Get the reference points. - const MatType& ReferencePoints() const { return referencePoints; } + const arma::Cube& ReferencePoints() const { return referencePoints; } //! Modify the reference points. - MatType& ReferencePoints() { return referencePoints; } + arma::Cube& ReferencePoints() { return referencePoints; } //! Retrieve value of lowerBound. const arma::vec& LowerBound() const { return lowerBound; } @@ -222,21 +222,20 @@ class NSGA3 * @param calculatedObjectives Vector to store calculated objectives. */ template, typename ...ArbitraryFunctionType> typename std::enable_if::type EvaluateObjectives(std::vector&, std::tuple&, - std::vector >&); + std::vector&); template, typename ...ArbitraryFunctionType> typename std::enable_if::type EvaluateObjectives(std::vector& population, std::tuple& objectives, - std::vector >& - calculatedObjectives); + std::vector& calculatedObjectives); /** * Reproduce candidates from the elite population to generate a new @@ -248,7 +247,6 @@ class NSGA3 * @param lowerBound Lower bound of the coordinates of the initial population. * @param upperBound Upper bound of the coordinates of the initial population. */ - template void BinaryTournamentSelection(std::vector& population, const MatType& lowerBound, const MatType& upperBound); @@ -262,7 +260,6 @@ class NSGA3 * @param parentA First parent from elite population. * @param parentB Second parent from elite population. */ - template void Crossover(MatType& childA, MatType& childB, const MatType& parentA, @@ -277,7 +274,6 @@ class NSGA3 * @param lowerBound Lower bound of the coordinates of the initial population. * @param upperBound Upper bound of the coordinates of the initial population. */ - template void Mutate(MatType& child, const MatType& lowerBound, const MatType& upperBound); @@ -292,11 +288,11 @@ class NSGA3 * @param ranks The assigned ranks, used for crowding distance based sorting. * @param calculatedObjectives The previously calculated objectives. */ - template + template void FastNonDominatedSort( std::vector >& fronts, std::vector& ranks, - std::vector >& calculatedObjectives); + std::vector& calculatedObjectives); /** * Normalizes the front given the extreme points in the current front. @@ -307,9 +303,9 @@ class NSGA3 * @param front The previously generated Pareto front. * @param extreme The indexes of the extreme points in the front. */ - void NormalizeFront( - std::vector >& calculatedObjectives, - arma::Col& normalization, + template + void NormalizeFront(std::vector& calculatedObjectives, + ColType& normalization, const std::vector& front, const arma::Row& extreme); @@ -320,9 +316,10 @@ class NSGA3 * @param calculatedObjectives The current population objectives. * @param front The front of the current generation. */ + template void FindExtremePoints(arma::Row& indexes, - std::vector >& calculatedObjectives, - const std::vector& front); + std::vector& calculatedObjectives, + const std::vector& front); /** * Finding the distance of each point in the front from the line formed @@ -333,26 +330,28 @@ class NSGA3 * @param front The front of the current generation(indices of population). * @param pointA The first point on the line. * @param pointB The second point on the line. - */ + */ + template void PointToLineDistance(arma::Row& distances, - std::vector >& calculatedObjectives, + std::vector& calculatedObjectives, const std::vector& front, - const arma::Col& pointA, - const arma::Col& pointB); + const ColType& pointA, + const ColType& pointB); /** * Finding the point in the reference front associated with the members in * the given front and also stores the distance between the two points. * - * @param index Vector containing the index of the point in the refrence directons associated. + * @param refIndex Vector containing the index of the point in the refrence directons associated. * @param dists Vector of distances from the corresponding point in the front to the associated reference direction. - * @param population The points of the currently generated population. - * @param front The index of points belonging to the given front. + * @param calculatedObjectives The points of the currently generated population. + * @param St The index of points belonging to the given front. */ - void Associate(arma::Row& refindex, + template + void Associate(arma::urowvec& refIndex, arma::Row& dists, - const std::vector& population, - const std::vector& front); + const std::vector& calculatedObjectives, + const std::vector& St); /** * Find the niche count for each reference direction. @@ -362,21 +361,25 @@ class NSGA3 * directons associated. */ void NicheCount(arma::Row& count, - const arma::Row& refIndex); + const arma::urowvec& refIndex); /** * The niche preserving operation to select the final points form the given front * aranges the front in descending order of priority for the top K points in the front. * * @param K The no. of remaining points to select from the given front for the next population. + * @param nicheCount The count of no. of points associated to each reference point in St. * @param refIndex The index of the rerference points associated with the in the given front. * @param dists The distances of th points in the front form their associated reference points line. * @param front The index of teh points within the population which are a part of the given front. + * @param population The set St (selected points). */ - void Niching(size_t K. - const arma::Row& refindex, + void Niching(size_t K, + arma::Row& nicheCount, + const arma::urowvec& refIndex, const arma::Row& dists, - const std::vector& front); + const std::vector& front, + std::vector& population); /** * Operator to check if one candidate Pareto-dominates the other. @@ -391,34 +394,12 @@ class NSGA3 * @param candidateQ The candidate being compared against. * @return true if candidateP Pareto dominates candidateQ, otherwise, false. */ - template + template bool Dominates( - std::vector >& calculatedObjectives, + std::vector& calculatedObjectives, size_t candidateP, size_t candidateQ); - /** - * The operator used in the crowding distance based sorting. - * - * If a candidates has a lower rank then it is preferred. - * Otherwise, if the ranks are equal then the candidate with the larger - * crowding distance is preferred. - * - * @param idxP The index of the first cadidate from the elite population being - * sorted. - * @param idxQ The index of the second cadidate from the elite population - * being sorted. - * @param ranks The previously calculated ranks. - * @param crowdingDistance The crowding distance for each individual in - * the population. - * @return true if the first candidate is preferred, otherwise, false. - */ - template - bool CrowdingOperator(size_t idxP, - size_t idxQ, - const std::vector& ranks, - const std::vector& crowdingDistance); - //! The number of objectives being optimised for. size_t numObjectives; @@ -458,7 +439,7 @@ class NSGA3 arma::cube paretoFront; //! The reference points. - MatType referencePoints; + arma::Cube referencePoints; //! A different representation of the Pareto front, for reverse compatibility //! purposes. This can be removed when ensmallen 3.x is released! (Along diff --git a/include/ensmallen_bits/nsga3/nsga3_impl.hpp b/include/ensmallen_bits/nsga3/nsga3_impl.hpp index 7cb04fdfc..c73fc756a 100644 --- a/include/ensmallen_bits/nsga3/nsga3_impl.hpp +++ b/include/ensmallen_bits/nsga3/nsga3_impl.hpp @@ -29,17 +29,17 @@ inline NSGA3::NSGA3(const MatType& referencePoints, const double epsilon = 1e-6, const arma::vec& lowerBound = arma::zeros(1, 1), const arma::vec& upperBound = arma::ones(1, 1)): - referencePoints(referencePoints), - numObjectives(0), - numVariables(0), - populationSize(populationSize), - maxGenerations(maxGenerations), - crossoverProb(crossoverProb), - distributionIndex(distributionIndex), - epsilon(epsilon), - eta(eta), - lowerBound(lowerBound), - upperBound(upperBound) + referencePoints(referencePoints), + numObjectives(0), + numVariables(0), + populationSize(populationSize), + maxGenerations(maxGenerations), + crossoverProb(crossoverProb), + distributionIndex(distributionIndex), + epsilon(epsilon), + eta(eta), + lowerBound(lowerBound), + upperBound(upperBound) { /* Nothing to do here. */ } template @@ -52,17 +52,17 @@ inline NSGA3::NSGA3(const MatType& referencePoints, const double epsilon = 1e-6, const double lowerBound = 0, const double upperBound = 1): - referencePoints(referencePoints), - numObjectives(0), - numVariables(0), - populationSize(populationSize), - maxGenerations(maxGenerations), - crossoverProb(crossoverProb), - distributionIndex(distributionIndex), - epsilon(epsilon), - eta(eta), - lowerBound(lowerBound * arma::ones(1, 1)), - upperBound(upperBound * arma::ones(1, 1)) + referencePoints(referencePoints), + numObjectives(0), + numVariables(0), + populationSize(populationSize), + maxGenerations(maxGenerations), + crossoverProb(crossoverProb), + distributionIndex(distributionIndex), + epsilon(epsilon), + eta(eta), + lowerBound(lowerBound * arma::ones(1, 1)), + upperBound(upperBound * arma::ones(1, 1)) { /* Nothing to do here. */ } //! Optimize the function. @@ -115,13 +115,13 @@ typename MatType::elem_type NSGA3::Optimize( // Population size reserved to 2 * populationSize + 1 to accommodate // for the size of intermediate candidate population. std::vector population; + std::vector tempPopulation; population.reserve(2 * populationSize + 1); + tempPopulation.reserve(populationSize); // Pareto fronts, initialized during non-dominated sorting. // Stores indices of population belonging to a certain front. std::vector > fronts; - // Initialised in SurvivalScoreAssignment. - std::vector survivalScore; // Initialised during non-dominated sorting. std::vector ranks; @@ -164,42 +164,56 @@ typename MatType::elem_type NSGA3::Optimize( ranks.resize(population.size()); FastNonDominatedSort(fronts, ranks, calculatedObjectives); + // S_t and P_t+1 declared. std::vector selectedPoints(fronts[0].begin(), fronts[0].end()); - while (selectedPoints.size() < populationSize) + std::vector nextPopulation(fronts[0].begin(), fronts[0].end()); + + size_t index = 0; + while (nextPopulation.size() + fronts.size() < populationSize) { - selectedPoints.insert(selectedPoints.end(), fronts[index].begin(), fronts[index].end()); + selectedPoints.insert(selectedPoints.end(), fronts[index].begin(), fronts[index].end()); + nextPopulation.insert(nextPopulation.end(), fronts[index].begin(), fronts[index].end()); + index++; } + selectedPoints.insert(selectedPoints.end(), fronts[++index].begin(), fronts[++index].end()); + arma::Col idealPoint(calculatedObjectives[selectedPoints[0]]); - for (size_t index = 1; index < selectedPoints.size(); index++) + for (index = 1; index < selectedPoints.size(); index++) { idealPoint = arma::min(idealPoint, calculatedObjectives[selectedPoints[index]]); } - arma::Col normalize(numObjectives, arma::fill::zeros); - - - // Sort based on survival score. - std::sort(population.begin(), population.end(), - [this, ranks, survivalScore, population] - (BaseMatType candidateP, BaseMatType candidateQ) - { - size_t idxP{}, idxQ{}; - for (size_t i = 0; i < population.size(); i++) - { - if (arma::approx_equal(population[i], candidateP, "absdiff", epsilon)) - idxP = i; - - if (arma::approx_equal(population[i], candidateQ, "absdiff", epsilon)) - idxQ = i; - } - - return SurvivalScoreOperator(idxP, idxQ, ranks, survivalScore); - } - ); + for (index = 0; index < selectedPoints.size(); index++) + { + calculatedObjectives[front[index]] = calculatedObjectives[front[index]] + - idealPoint; + } - // Yield a new population P_{t+1} of size populationSize. - // Discards unfit population from the R_{t} to yield P_{t+1}. - population.resize(populationSize); + arma::Col normalize(numObjectives, arma::fill::zeros); + arma::Row extreme(numObjectives, arma::fill::zeros); + + // Find the extreme points and normalize S_t + FindExtremePoints>(extreme, calculatedObjectives, + selectedPoints); + NormalizeFront>(calculatedObjectives, normalize, + selectedPoints, extreme); + + // Find the associated reference directions to the selected points. + arma::urowvec refIndex(selectedPoints.size()); + arma::Row dists(selectedPoints.size()); + Associate>(refIndex, dists, calculatedObjectives, + selectedPoints); + + // Calculate the niche count of S_t and performing the niching operation. + arma::Row count(referencePoints.n_slices); + NicheCount(count, refIndex); + Niching(K, nicheCount, refIndex, dists, front, nextPopulation) + + for (size_t i : nextPopulation) + { + tempPopulation.push_back(population[i]); + } + population = tempPopulation; terminate |= Callback::GenerationalStepTaken(*this, objectives, iterate, calculatedObjectives, fronts, callbacks...); @@ -243,27 +257,28 @@ typename MatType::elem_type NSGA3::Optimize( } //! No objectives to evaluate. +template template typename std::enable_if::type -AGEMOEA::EvaluateObjectives( +NSGA3::EvaluateObjectives( std::vector&, std::tuple&, - std::vector >&) + std::vector >&) { // Nothing to do here. } //! Evaluate the objectives for the entire population. template typename std::enable_if::type -AGEMOEA::EvaluateObjectives( +NSGA3::EvaluateObjectives( std::vector& population, std::tuple& objectives, - std::vector >& calculatedObjectives) + std::vector& calculatedObjectives) { for (size_t i = 0; i < populationSize; i++) { @@ -275,9 +290,10 @@ AGEMOEA::EvaluateObjectives( //! Reproduce and generate new candidates. template -inline void AGEMOEA::BinaryTournamentSelection(std::vector& population, - const MatType& lowerBound, - const MatType& upperBound) +inline void NSGA3::BinaryTournamentSelection( + std::vector& population, + const MatType& lowerBound, + const MatType& upperBound) { std::vector children; @@ -318,13 +334,13 @@ inline void AGEMOEA::BinaryTournamentSelection(std::vector& population, } //! Perform simulated binary crossover (SBX) of genes for the children. -template -inline void AGEMOEA::Crossover(MatType& childA, - MatType& childB, - const MatType& parentA, - const MatType& parentB, - const MatType& lowerBound, - const MatType& upperBound) +template +inline void NSGA3::Crossover(MatType& childA, + MatType& childB, + const MatType& parentA, + const MatType& parentB, + const MatType& lowerBound, + const MatType& upperBound) { //! Generates a child from two parent individuals // according to the polynomial probability distribution. @@ -377,10 +393,10 @@ inline void AGEMOEA::Crossover(MatType& childA, //! Perform Polynomial mutation of the candidate. template -inline void AGEMOEA::Mutate(MatType& candidate, - double mutationRate, - const MatType& lowerBound, - const MatType& upperBound) +inline void NSGA3::Mutate(MatType& candidate, + double mutationRate, + const MatType& lowerBound, + const MatType& upperBound) { const size_t numVariables = candidate.n_rows; for (size_t geneIdx = 0; geneIdx < numVariables; ++geneIdx) @@ -416,11 +432,12 @@ inline void AGEMOEA::Mutate(MatType& candidate, } template -inline void AGEMOEA::NormalizeFront( - std::vector >& calculatedObjectives, - arma::Col& normalization, - const std::vector& front, - const arma::Row& extreme) +template +inline void NSGA3::NormalizeFront( + std::vector& calculatedObjectives, + ColType& normalization, + const std::vector& front, + const arma::Row& extreme) { arma::Mat vectorizedObjectives(numObjectives, front.size()); for (size_t i = 0; i < front.size(); i++) @@ -433,15 +450,15 @@ inline void AGEMOEA::NormalizeFront( normalization = arma::max(vectorizedObjectives, 1); return; } - arma::Col temp; + ColType temp; arma::uvec unique = arma::find_unique(extreme); if (extreme.n_elem != unique.n_elem) { normalization = arma::max(vectorizedObjectives, 1); return; } - arma::Col one(front.size(), arma::fill::ones); - arma::Col hyperplane = arma::solve( + ColType one(front.size(), arma::fill::ones); + ColType hyperplane = arma::solve( vectorizedObjectives.t(), one); if (hyperplane.has_inf() || hyperplane.has_nan() || (arma::accu(hyperplane < 0.0) > 0)) { @@ -460,9 +477,10 @@ inline void AGEMOEA::NormalizeFront( //! Find the index of the of the extreme points in the given front. template +template void NSGA::FindExtremePoints(arma::Row& indexes, - std::vector >& calculatedObjectives, - const std::vector& front) + std::vector& calculatedObjectives, + const std::vector& front) { typedef typename MatType::elem_type ElemType; @@ -479,7 +497,7 @@ void NSGA::FindExtremePoints(arma::Row& indexes, arma::Row dists; for (size_t i = 0; i < numObjectives; i++) { - PointToLineDistance(dists, calculatedObjectives, front, z, W.col(i)); + PointToLineDistance(dists, calculatedObjectives, front, z, W.col(i)); for (size_t j = 0; j < front.size(); j++) if (selected[j]){dists[j] = arma::datum::inf;} indexes[i] = dists.index_min(); @@ -489,16 +507,18 @@ void NSGA::FindExtremePoints(arma::Row& indexes, //! Find the distance of a front from a line formed by two points. template -void NSGA3::PointToLineDistance(arma::Row& distances, - std::vector >& calculatedObjectives, - const std::vector& front, - const arma::Col& pointA, - const arma::Col& pointB) +template +void NSGA3::PointToLineDistance( + arma::Row& distances, + std::vector& calculatedObjectives, + const std::vector& front, + const ColType& pointA, + const ColType& pointB) { typedef typename MatType::elem_type ElemType; arma::Row distancesTemp(front.size()); - arma::Col ba = pointB - pointA; - arma::Col pa; + ColType ba = pointB - pointA; + ColType pa; for (size_t i = 0; i < front.size(); i++) { @@ -513,10 +533,11 @@ void NSGA3::PointToLineDistance(arma::Row& //! Sort population into Pareto fronts. template -inline void AGEMOEA::FastNonDominatedSort( +template +inline void NSGA3::FastNonDominatedSort( std::vector >& fronts, std::vector& ranks, - std::vector >& calculatedObjectives) + std::vector& calculatedObjectives) { std::map dominationCount; std::map > dominated; @@ -572,10 +593,84 @@ inline void AGEMOEA::FastNonDominatedSort( fronts.pop_back(); } +template +void NSGA3::Niching(size_t K, + arma::Row& nicheCount, + const arma::urowvec& refIndex, + const arma::Row& dists, + const std::vector& front, + std::vector& population) // population excluding current front. +{ + typedef typename MatType::elem_type elemType; + arma::Row popMask(population.n_elem); + size_t k = 0; + while (k < K) + { + size_t jMin = arma::index_min(nicheCount); + std::vector I; + for (size_t i : front) + { + if(refIndex[i] == jMin && popMask[i]) + I.push_back(i); + } + if (I.size() != 0) + { + size_t min = 0; + if(nicheCount[jMin] == 0) + { + for (size_t i = 0; i < I.size(); i++) + { + if(dists[I[i]] < dists[I[min]]) + { + min = i; + } + } + } + population.append(I[min]); + + nicheCount[jMin] += 1; + frontMask[I[i]] = false; + k += 1; + } + else + { + nicheCount[jMin] = arma::datum::inf; + } + } +} + + +template +template +void NSGA3::Associate(arma::urowvec& refindex, + arma::Row& dists, + const std::vector& calculatedObjectives, + const std::vector& St) +{ + arma::Mat d(refindex.n_elem, St.size()); + for (size_t i = 0; i < referencePoints.n_slices; i++) + { + PointToLineDistance(d.row(i), calculatedObjectives, St, , W.col(i)); + } + refIndex = arma::index_min(d, 0); + dists = arma::min(d, 0); +} + +template +void NSGA3::NicheCount(arma::Row& count, + const arma::urowvec& refIndex) +{ + for (size_t i = 0; i < refIndex.n_elem; i++) + { + count[refIndex[i]] += 1; + } +} + //! Check if a candidate Pareto dominates another candidate. template -inline bool AGEMOEA::Dominates( - std::vector >& calculatedObjectives, +template +inline bool NSGA3::Dominates( + std::vector& calculatedObjectives, size_t candidateP, size_t candidateQ) { From 8dba616fae7a10f222e5ca75383332975792c9ef Mon Sep 17 00:00:00 2001 From: IWNMWE Date: Sat, 3 Aug 2024 01:30:31 +0530 Subject: [PATCH 03/11] nsga3 definition --- include/ensmallen.hpp | 1 + include/ensmallen_bits/nsga3/nsga3.hpp | 52 +++++++++------ include/ensmallen_bits/nsga3/nsga3_impl.hpp | 74 +++++++++++---------- 3 files changed, 72 insertions(+), 55 deletions(-) diff --git a/include/ensmallen.hpp b/include/ensmallen.hpp index a6ca139dd..f596c7734 100644 --- a/include/ensmallen.hpp +++ b/include/ensmallen.hpp @@ -120,6 +120,7 @@ #include "ensmallen_bits/agemoea/agemoea.hpp" #include "ensmallen_bits/moead/moead.hpp" #include "ensmallen_bits/nsga2/nsga2.hpp" +#include "ensmallen_bits/nsga3/nsga3.hpp" #include "ensmallen_bits/padam/padam.hpp" #include "ensmallen_bits/parallel_sgd/parallel_sgd.hpp" #include "ensmallen_bits/pso/pso.hpp" diff --git a/include/ensmallen_bits/nsga3/nsga3.hpp b/include/ensmallen_bits/nsga3/nsga3.hpp index 0831e9257..c95f3389a 100644 --- a/include/ensmallen_bits/nsga3/nsga3.hpp +++ b/include/ensmallen_bits/nsga3/nsga3.hpp @@ -79,13 +79,13 @@ class NSGA3 * @param lowerBound Lower bound of the coordinates of the initial population. * @param upperBound Upper bound of the coordinates of the initial population. */ - NSGA3(const arma::Cube& referencePoints, + NSGA3(const arma::Cube& referencePoints, const size_t populationSize = 100, const size_t maxGenerations = 2000, const double crossoverProb = 0.6, - const double mutationProb = 0.3, - const double mutationStrength = 1e-3, + const double distributionIndex = 20, const double epsilon = 1e-6, + const double eta = 20, const arma::vec& lowerBound = arma::zeros(1, 1), const arma::vec& upperBound = arma::ones(1, 1)); @@ -109,13 +109,13 @@ class NSGA3 * @param lowerBound Lower bound of the coordinates of the initial population. * @param upperBound Upper bound of the coordinates of the initial population. */ - NSGA3(const arma::Cube& referencePoints, + NSGA3(const arma::Cube& referencePoints, const size_t populationSize = 100, const size_t maxGenerations = 2000, const double crossoverProb = 0.6, - const double mutationProb = 0.3, - const double mutationStrength = 1e-3, + const double distributionIndex = 20, const double epsilon = 1e-6, + const double eta = 20, const double lowerBound = 0, const double upperBound = 1); @@ -153,26 +153,21 @@ class NSGA3 //! Modify the crossover rate. double& CrossoverRate() { return crossoverProb; } - //! Get the mutation probability. - double MutationProbability() const { return mutationProb; } - //! Modify the mutation probability. - double& MutationProbability() { return mutationProb; } + //! Retrieve value of the distribution index. + double DistributionIndex() const { return distributionIndex; } + //! Modify the value of the distribution index. + double& DistributionIndex() { return distributionIndex; } - //! Get the mutation strength. - double MutationStrength() const { return mutationStrength; } - //! Modify the mutation strength. - double& MutationStrength() { return mutationStrength; } + //! Retrieve value of eta. + double Eta() const { return eta; } + //! Modify the value of eta. + double& Eta() { return eta; } //! Get the tolerance. double Epsilon() const { return epsilon; } //! Modify the tolerance. double& Epsilon() { return epsilon; } - //! Get the reference points. - const arma::Cube& ReferencePoints() const { return referencePoints; } - //! Modify the reference points. - arma::Cube& ReferencePoints() { return referencePoints; } - //! Retrieve value of lowerBound. const arma::vec& LowerBound() const { return lowerBound; } //! Modify value of lowerBound. @@ -191,6 +186,11 @@ class NSGA3 //! `Optimize()` has been called. const arma::cube& ParetoFront() const { return paretoFront; } + //! Get the reference points. + const arma::Cube& ReferencePoints() const { return referencePoints; } + //! Modify the reference points. + arma::Cube& ReferencePoints() { return referencePoints; } + /** * Retrieve the best front (the Pareto frontier). This returns an empty * vector until `Optimize()` has been called. Note that this function is @@ -263,7 +263,9 @@ class NSGA3 void Crossover(MatType& childA, MatType& childB, const MatType& parentA, - const MatType& parentB); + const MatType& parentB, + const MatType& lowerBound, + const MatType& upperBound); /** * Mutate the coordinates for a candidate. @@ -274,7 +276,8 @@ class NSGA3 * @param lowerBound Lower bound of the coordinates of the initial population. * @param upperBound Upper bound of the coordinates of the initial population. */ - void Mutate(MatType& child, + void Mutate(MatType& candidate, + double mutationRate, const MatType& lowerBound, const MatType& upperBound); @@ -421,9 +424,16 @@ class NSGA3 //! Strength of the mutation. double mutationStrength; + //! The crowding degree of the mutation. Higher value produces a mutant + //! resembling its parent. + double distributionIndex; + //! The tolerance for termination. double epsilon; + //! The distance parameters of the crossover distribution. + double eta; + //! Lower bound of the initial swarm. arma::vec lowerBound; diff --git a/include/ensmallen_bits/nsga3/nsga3_impl.hpp b/include/ensmallen_bits/nsga3/nsga3_impl.hpp index c73fc756a..32e533ed9 100644 --- a/include/ensmallen_bits/nsga3/nsga3_impl.hpp +++ b/include/ensmallen_bits/nsga3/nsga3_impl.hpp @@ -20,15 +20,16 @@ namespace ens { template -inline NSGA3::NSGA3(const MatType& referencePoints, - const size_t populationSize = 100, - const size_t maxGenerations = 2000, - const double crossoverProb = 0.6, - const double mutationProb = 0.3, - const double mutationStrength = 1e-3, - const double epsilon = 1e-6, - const arma::vec& lowerBound = arma::zeros(1, 1), - const arma::vec& upperBound = arma::ones(1, 1)): +inline NSGA3::NSGA3( + const arma::Cube& referencePoints, + const size_t populationSize, + const size_t maxGenerations, + const double crossoverProb, + const double distributionIndex, + const double epsilon, + const double eta, + const arma::vec& lowerBound, + const arma::vec& upperBound): referencePoints(referencePoints), numObjectives(0), numVariables(0), @@ -43,15 +44,16 @@ inline NSGA3::NSGA3(const MatType& referencePoints, { /* Nothing to do here. */ } template -inline NSGA3::NSGA3(const MatType& referencePoints, - const size_t populationSize = 100, - const size_t maxGenerations = 2000, - const double crossoverProb = 0.6, - const double mutationProb = 0.3, - const double mutationStrength = 1e-3, - const double epsilon = 1e-6, - const double lowerBound = 0, - const double upperBound = 1): +inline NSGA3::NSGA3( + const arma::Cube& referencePoints, + const size_t populationSize, + const size_t maxGenerations, + const double crossoverProb, + const double distributionIndex, + const double epsilon, + const double eta, + const double lowerBound, + const double upperBound): referencePoints(referencePoints), numObjectives(0), numVariables(0), @@ -176,6 +178,7 @@ typename MatType::elem_type NSGA3::Optimize( index++; } selectedPoints.insert(selectedPoints.end(), fronts[++index].begin(), fronts[++index].end()); + size_t lastFront = index; arma::Col idealPoint(calculatedObjectives[selectedPoints[0]]); for (index = 1; index < selectedPoints.size(); index++) @@ -185,8 +188,8 @@ typename MatType::elem_type NSGA3::Optimize( for (index = 0; index < selectedPoints.size(); index++) { - calculatedObjectives[front[index]] = calculatedObjectives[front[index]] - - idealPoint; + calculatedObjectives[selectedPoints[index]] = + calculatedObjectives[selectedPoints[index]] - idealPoint; } arma::Col normalize(numObjectives, arma::fill::zeros); @@ -207,7 +210,8 @@ typename MatType::elem_type NSGA3::Optimize( // Calculate the niche count of S_t and performing the niching operation. arma::Row count(referencePoints.n_slices); NicheCount(count, refIndex); - Niching(K, nicheCount, refIndex, dists, front, nextPopulation) + Niching(selectedPoints.size() - nextPopulation.size(), count, refIndex, + dists, fronts[lastFront], nextPopulation); for (size_t i : nextPopulation) { @@ -265,12 +269,13 @@ typename std::enable_if::type NSGA3::EvaluateObjectives( std::vector&, std::tuple&, - std::vector >&) + std::vector&) { // Nothing to do here. } //! Evaluate the objectives for the entire population. +template template @@ -394,9 +399,9 @@ inline void NSGA3::Crossover(MatType& childA, //! Perform Polynomial mutation of the candidate. template inline void NSGA3::Mutate(MatType& candidate, - double mutationRate, - const MatType& lowerBound, - const MatType& upperBound) + double mutationRate, + const MatType& lowerBound, + const MatType& upperBound) { const size_t numVariables = candidate.n_rows; for (size_t geneIdx = 0; geneIdx < numVariables; ++geneIdx) @@ -478,7 +483,7 @@ inline void NSGA3::NormalizeFront( //! Find the index of the of the extreme points in the given front. template template -void NSGA::FindExtremePoints(arma::Row& indexes, +void NSGA3::FindExtremePoints(arma::Row& indexes, std::vector& calculatedObjectives, const std::vector& front) { @@ -599,10 +604,9 @@ void NSGA3::Niching(size_t K, const arma::urowvec& refIndex, const arma::Row& dists, const std::vector& front, - std::vector& population) // population excluding current front. + std::vector& population) { - typedef typename MatType::elem_type elemType; - arma::Row popMask(population.n_elem); + arma::Row popMask(population.size()); size_t k = 0; while (k < K) { @@ -626,10 +630,10 @@ void NSGA3::Niching(size_t K, } } } - population.append(I[min]); + population.push_back(I[min]); nicheCount[jMin] += 1; - frontMask[I[i]] = false; + popMask[I[min]] = false; k += 1; } else @@ -642,15 +646,17 @@ void NSGA3::Niching(size_t K, template template -void NSGA3::Associate(arma::urowvec& refindex, +void NSGA3::Associate(arma::urowvec& refIndex, arma::Row& dists, const std::vector& calculatedObjectives, const std::vector& St) { - arma::Mat d(refindex.n_elem, St.size()); + arma::Mat d(refIndex.n_elem, St.size()); + ColType zero(arma::size(calculatedObjectives),arma::fill::zeros); for (size_t i = 0; i < referencePoints.n_slices; i++) { - PointToLineDistance(d.row(i), calculatedObjectives, St, , W.col(i)); + PointToLineDistance(d.row(i), calculatedObjectives, St, + referencePoints.slice(i), zero); } refIndex = arma::index_min(d, 0); dists = arma::min(d, 0); From 279fcf971417f9c685c6beafb53b479b75d8a002 Mon Sep 17 00:00:00 2001 From: IWNMWE Date: Thu, 22 Aug 2024 02:04:10 +0530 Subject: [PATCH 04/11] Add NSGA3 normalization --- include/ensmallen.hpp | 1 + .../ensmallen_bits/nsga3/normalization.hpp | 235 ++++++++++++++++++ 2 files changed, 236 insertions(+) create mode 100644 include/ensmallen_bits/nsga3/normalization.hpp diff --git a/include/ensmallen.hpp b/include/ensmallen.hpp index f596c7734..8963d15c0 100644 --- a/include/ensmallen.hpp +++ b/include/ensmallen.hpp @@ -120,6 +120,7 @@ #include "ensmallen_bits/agemoea/agemoea.hpp" #include "ensmallen_bits/moead/moead.hpp" #include "ensmallen_bits/nsga2/nsga2.hpp" +#include "ensmallen_bits/nsga3/normalization.hpp" #include "ensmallen_bits/nsga3/nsga3.hpp" #include "ensmallen_bits/padam/padam.hpp" #include "ensmallen_bits/parallel_sgd/parallel_sgd.hpp" diff --git a/include/ensmallen_bits/nsga3/normalization.hpp b/include/ensmallen_bits/nsga3/normalization.hpp new file mode 100644 index 000000000..01f0bd436 --- /dev/null +++ b/include/ensmallen_bits/nsga3/normalization.hpp @@ -0,0 +1,235 @@ +/** + * @file normalization.hpp + * @author Satyam Shukla + * + * The Optimized normalization technique as described in Investigating the + * normalization procedure of nsga-iii. + * + * ensmallen is free software; you may redistribute it and/or modify it under + * the terms of the 3-clause BSD license. You should have received a copy of + * the 3-clause BSD license along with ensmallen. If not, see + * http://www.opensource.org/licenses/BSD-3-Clause for more information. + */ +#ifndef ENSMALLEN_NSGA3_NORMALIZATION_HPP +#define ENSMALLEN_NSGA3_NORMALIZATION_HPP + +namespace ens { + +/** + * + * This normalization technique is an improved version of the algorithm mentioned + * in the NSGA III paper. + * + * This class solves the problem of negative intercepts and non unique hyperplane + * by using worst point estimation as a backup value when ther are any abnormal values + * in the nadir point estimation. + * + * For more information, see the following: + * + * @code + * @inproceedings{blank2019investigating, + * title={Investigating the normalization procedure of NSGA-III}, + * author={Blank, Julian and Deb, Kalyanmoy and Roy, Proteek Chandan}, + * booktitle={International Conference on Evolutionary Multi-Criterion Optimization}, + * pages={229--240}, + * year={2019}, + * organization={Springer} + * } + * @endcode + */ +template +class Normalization +{ + public: + + typedef typename MatType::elem_type ElemType; + + /** + * @param dimension The no. of elements in a single point. + */ + Normalization(size_t dimensions = 3): + dimensions(dimensions), + idealPoint(arma::Col(dimensions, + arma::fill::value(arma::datum::inf))), + worstPoint(arma::Col(dimensions, + arma::fill::value(-1 * arma::datum::inf))) + {/* Nothing to do here */} + + /** + * @param calculatedObjectives The given population. + * @param indexes The index of the lements of the front + */ + void update(const std::vector>& calculatedObjectives, + const std::vector indexes) + { + // Calculate the front worst, the population worst and the ideal point. + arma::Col worstOfFront(dimensions, + arma::fill::value(-1 * arma::datum::inf)); + arma::Col worstOfPop(dimensions, + arma::fill::value(-1 * arma::datum::inf)); + + for (const arma::Col member: calculatedObjectives) + { + idealPoint = arma::min(idealPoint, member); + worstPoint = arma::max(worstPoint, member); + worstOfPop = arma::max(worstOfPop, member); + } + + for (size_t index : indexes) + { + worstOfFront = arma::max(worstOfFront, calculatedObjectives[index]); + } + + arma::Mat f(dimensions, dimensions); + + // Find the extremes. + FindExtremes(calculatedObjectives, indexes, f); + vectorizedExtremes = f; + + // Update the nadir point. + GetNadirPoint(calculatedObjectives, indexes, worstOfPop, worstOfFront); + + } + + /** + * @param calculatedObjectives The given population. + * @param indexes The index of the lements of the front + * @param f The matrix to store the vector extremes. + * @param useCurrentExtremes If the previously calculaeted extremes should + * be considered. + */ + void FindExtremes(const std::vector>& calculatedObjectives, + const std::vector& indexes, + arma::Mat& f, + bool useCurrentExtremes = true) + { + arma::Mat W(dimensions, dimensions, arma::fill::eye); + W = W + (W == 0) * 1e6; + arma::Mat vectorizedObjectives(dimensions, indexes.size()); + + for (size_t i = 0; i < indexes.size(); i++) + { + vectorizedObjectives.col(i) = calculatedObjectives[indexes[i]]; + } + + if (useCurrentExtremes) + { + vectorizedObjectives = arma::join_rows(vectorizedObjectives, + vectorizedExtremes); + } + vectorizedObjectives.each_col() -= idealPoint; + vectorizedObjectives = vectorizedObjectives + (vectorizedObjectives < 1e-3) * 0.; + + //Calculate ASF score and get the extreme vectors. + arma::Mat asfScore(vectorizedObjectives.n_cols, dimensions); + for(size_t i = 0; i < vectorizedObjectives.n_cols; i++) + { + for(size_t j = 0; j < dimensions; j++) + { + asfScore(i, j) = arma::max(W.col(j) % vectorizedObjectives.col(i)); + } + } + arma::urowvec extremes = arma::index_min(asfScore); + for(size_t i = 0; i < dimensions; i++) + { + if(extremes(i) >= indexes.size()) + { + f.col(i) = vectorizedExtremes.col(extremes(i) - indexes.size()); + } + else + { + f.col(i) = calculatedObjectives[indexes[extremes(i)]]; + } + } + } + + /** + * @param calculatedObjectives The given population. + * @param indexes The index of the lements of the front + * @param worstOfFront Worst point of the given front. + * @param worstOfPop Worst point of the population. + */ + void GetNadirPoint(const std::vector>& calculatedObjectives, + const std::vector& indexes, + const arma::Col worstOfFront, + const arma::Col worstOfPop) + { + try + { + arma::Mat M = vectorizedExtremes; + M.each_col() -= idealPoint; + arma::Col b(dimensions, arma::fill::ones); + arma::Col hyperplane = arma::solve(M.t(), b); + + if (hyperplane.has_inf() || hyperplane.has_nan() || (arma::accu(hyperplane < 0.0) > 0)) + { + throw 1024; + } + + arma::Col intercepts = 1.0 / hyperplane; + + if(arma::accu(arma::abs((M.t() * hyperplane) - b) > 1e-8) || + arma::accu(intercepts < 1e-6) > 0 || intercepts.has_inf() || + intercepts.has_nan()) + { + throw 1025; + } + + nadirPoint = idealPoint + intercepts; + nadirPoint = nadirPoint % (nadirPoint <= worstPoint); + nadirPoint += (nadirPoint == 0) % worstPoint; + } + catch(...) + { + nadirPoint = worstOfFront; + } + nadirPoint = ((nadirPoint - idealPoint) >= 1e-6) % nadirPoint + + ((nadirPoint - idealPoint) < 1e-6) % worstOfPop; + } + + //! Retrieve value of dimensions. + size_t Dimensions() const { return dimensions; } + //! Modify value of dimensions. + size_t& Dimensions() {return dimensions;} + + //! Retrieve value of ideal point. + arma::Col IdealPoint() const {return idealPoint; } + //! Modify value of ideal point. + arma::Col& IdealPoint() {return idealPoint; } + + //! Retrieve value of worst point. + arma::Col WorstPoint() const {return worstPoint; } + //! Modify value of worst point. + arma::Col& WorstPoint() {return worstPoint; } + + //! Retrieve value of nadir point. + arma::Col NadirPoint() const {return nadirPoint; } + //! Modify value of nadir point. + arma::Col& NadirPoint() {return nadirPoint; } + + //! Retrieve value of extreme points. + arma::Mat VectorizedExtremes() const {return vectorizedExtremes; } + //! Modify value of extreme points. + arma::Mat& VectorizedExtremes() {return vectorizedExtremes; } + + private: + + size_t dimensions; + + //The ideal point of the current population . + //! includes previous ideal point in calculations as well. + arma::Col idealPoint; + + // The worst point in the current population. + //! includes previous worst point in calculations as well. + arma::Col worstPoint; + + // The nadir point of the previous front. + arma::Col nadirPoint; + + // A Matrix containing the extreme vectors as columns. + arma::Mat vectorizedExtremes; +}; +} + +#endif From 23ef8af8a45b61ff63541ad8af0a0a543123e581 Mon Sep 17 00:00:00 2001 From: IWNMWE Date: Thu, 22 Aug 2024 02:04:32 +0530 Subject: [PATCH 05/11] Integrate normalization in nsga3 --- include/ensmallen_bits/nsga3/nsga3.hpp | 98 ++--- include/ensmallen_bits/nsga3/nsga3_impl.hpp | 348 ++++++++---------- .../problems/dtlz/dtlz1_function.hpp | 4 - 3 files changed, 186 insertions(+), 264 deletions(-) diff --git a/include/ensmallen_bits/nsga3/nsga3.hpp b/include/ensmallen_bits/nsga3/nsga3.hpp index c95f3389a..d37fe1715 100644 --- a/include/ensmallen_bits/nsga3/nsga3.hpp +++ b/include/ensmallen_bits/nsga3/nsga3.hpp @@ -56,7 +56,7 @@ namespace ens { * see the documentation on function types included with this distribution or * on the ensmallen website. */ -template +template class NSGA3 { public: @@ -74,17 +74,14 @@ class NSGA3 * @param crossoverProb The probability that a crossover will occur. * @param mutationProb The probability that a mutation will occur. * @param mutationStrength The strength of the mutation. - * @param epsilon The minimum difference required to distinguish between - * candidate solutions. * @param lowerBound Lower bound of the coordinates of the initial population. * @param upperBound Upper bound of the coordinates of the initial population. */ - NSGA3(const arma::Cube& referencePoints, + NSGA3(const arma::Mat& referencePoints, const size_t populationSize = 100, const size_t maxGenerations = 2000, const double crossoverProb = 0.6, const double distributionIndex = 20, - const double epsilon = 1e-6, const double eta = 20, const arma::vec& lowerBound = arma::zeros(1, 1), const arma::vec& upperBound = arma::ones(1, 1)); @@ -104,17 +101,14 @@ class NSGA3 * @param crossoverProb The probability that a crossover will occur. * @param mutationProb The probability that a mutation will occur. * @param mutationStrength The strength of the mutation. - * @param epsilon The minimum difference required to distinguish between - * candidate solutions. * @param lowerBound Lower bound of the coordinates of the initial population. * @param upperBound Upper bound of the coordinates of the initial population. */ - NSGA3(const arma::Cube& referencePoints, + NSGA3(const arma::Mat& referencePoints, const size_t populationSize = 100, const size_t maxGenerations = 2000, const double crossoverProb = 0.6, const double distributionIndex = 20, - const double epsilon = 1e-6, const double eta = 20, const double lowerBound = 0, const double upperBound = 1); @@ -124,6 +118,7 @@ class NSGA3 * starting point. The output is the best generated front. * * @tparam ArbitraryFunctionType std::tuple of multiple objectives. + * @tparam MatType Type of matrix to optimize. * @tparam CallbackTypes Types of callback functions. * @param objectives Vector of objective functions to optimize for. * @param iterate Starting point. @@ -131,12 +126,13 @@ class NSGA3 * @return MatType::elem_type The minimum of the accumulated sum over the * objective values in the best front. */ - template - typename MatType::elem_type Optimize( - std::tuple& objectives, - MatType& iterate, - CallbackTypes&&... callbacks); + typename MatType::elem_type Optimize( + std::tuple& objectives, + MatType& iterate, + CallbackTypes&&... callbacks); //! Get the population size. size_t PopulationSize() const { return populationSize; } @@ -163,11 +159,6 @@ class NSGA3 //! Modify the value of eta. double& Eta() { return eta; } - //! Get the tolerance. - double Epsilon() const { return epsilon; } - //! Modify the tolerance. - double& Epsilon() { return epsilon; } - //! Retrieve value of lowerBound. const arma::vec& LowerBound() const { return lowerBound; } //! Modify value of lowerBound. @@ -187,9 +178,9 @@ class NSGA3 const arma::cube& ParetoFront() const { return paretoFront; } //! Get the reference points. - const arma::Cube& ReferencePoints() const { return referencePoints; } + const arma::Mat& ReferencePoints() const { return referencePoints; } //! Modify the reference points. - arma::Cube& ReferencePoints() { return referencePoints; } + arma::Mat& ReferencePoints() { return referencePoints; } /** * Retrieve the best front (the Pareto frontier). This returns an empty @@ -211,7 +202,6 @@ class NSGA3 return rcFront; } - private: /** * Evaluate objectives for the elite population. * @@ -222,31 +212,32 @@ class NSGA3 * @param calculatedObjectives Vector to store calculated objectives. */ template, + typename MatType, typename ...ArbitraryFunctionType> typename std::enable_if::type EvaluateObjectives(std::vector&, std::tuple&, - std::vector&); + std::vector >&); template, + typename MatType, typename ...ArbitraryFunctionType> typename std::enable_if::type EvaluateObjectives(std::vector& population, std::tuple& objectives, - std::vector& calculatedObjectives); + std::vector >& + calculatedObjectives); /** * Reproduce candidates from the elite population to generate a new * population. * - * @tparam MatType Type of matrix to optimize. + * @tparam BaseMatType Type of matrix to optimize. * @param population The elite population. - * @param objectives The set of objectives. * @param lowerBound Lower bound of the coordinates of the initial population. * @param upperBound Upper bound of the coordinates of the initial population. */ + template void BinaryTournamentSelection(std::vector& population, const MatType& lowerBound, const MatType& upperBound); @@ -259,7 +250,10 @@ class NSGA3 * @param childB Another newly generated candidate. * @param parentA First parent from elite population. * @param parentB Second parent from elite population. + * @param lowerBound Lower Bound of the offspring. + * @param upperBound Upper Boundn of the offspring. */ + template void Crossover(MatType& childA, MatType& childB, const MatType& parentA, @@ -271,11 +265,12 @@ class NSGA3 * Mutate the coordinates for a candidate. * * @tparam MatType Type of matrix to optimize. - * @param child The candidate whose coordinates are being modified. - * @param objectives The set of objectives. + * @param candidate The candidate whose coordinates are being modified. + * @param mutationRate The probability of mutation. * @param lowerBound Lower bound of the coordinates of the initial population. * @param upperBound Upper bound of the coordinates of the initial population. */ + template void Mutate(MatType& candidate, double mutationRate, const MatType& lowerBound, @@ -296,33 +291,6 @@ class NSGA3 std::vector >& fronts, std::vector& ranks, std::vector& calculatedObjectives); - - /** - * Normalizes the front given the extreme points in the current front. - * - * @tparam The type of population datapoints. - * @param calculatedObjectives The current population evaluated objectives. - * @param normalization The normalizing vector. - * @param front The previously generated Pareto front. - * @param extreme The indexes of the extreme points in the front. - */ - template - void NormalizeFront(std::vector& calculatedObjectives, - ColType& normalization, - const std::vector& front, - const arma::Row& extreme); - - /** - * Finding the indexes of the extreme points in the front. - * - * @param indexes vector containing the slected indexes. - * @param calculatedObjectives The current population objectives. - * @param front The front of the current generation. - */ - template - void FindExtremePoints(arma::Row& indexes, - std::vector& calculatedObjectives, - const std::vector& front); /** * Finding the distance of each point in the front from the line formed @@ -335,8 +303,8 @@ class NSGA3 * @param pointB The second point on the line. */ template - void PointToLineDistance(arma::Row& distances, - std::vector& calculatedObjectives, + void PointToLineDistance(arma::Row& distances, + const std::vector& calculatedObjectives, const std::vector& front, const ColType& pointA, const ColType& pointB); @@ -352,7 +320,7 @@ class NSGA3 */ template void Associate(arma::urowvec& refIndex, - arma::Row& dists, + arma::Row& dists, const std::vector& calculatedObjectives, const std::vector& St); @@ -364,7 +332,8 @@ class NSGA3 * directons associated. */ void NicheCount(arma::Row& count, - const arma::urowvec& refIndex); + const arma::urowvec& refIndex, + const std::vector& nextPopulation); /** * The niche preserving operation to select the final points form the given front @@ -380,7 +349,7 @@ class NSGA3 void Niching(size_t K, arma::Row& nicheCount, const arma::urowvec& refIndex, - const arma::Row& dists, + const arma::Row& dists, const std::vector& front, std::vector& population); @@ -428,9 +397,6 @@ class NSGA3 //! resembling its parent. double distributionIndex; - //! The tolerance for termination. - double epsilon; - //! The distance parameters of the crossover distribution. double eta; @@ -449,7 +415,7 @@ class NSGA3 arma::cube paretoFront; //! The reference points. - arma::Cube referencePoints; + arma::Mat referencePoints; //! A different representation of the Pareto front, for reverse compatibility //! purposes. This can be removed when ensmallen 3.x is released! (Along diff --git a/include/ensmallen_bits/nsga3/nsga3_impl.hpp b/include/ensmallen_bits/nsga3/nsga3_impl.hpp index 32e533ed9..ebe513509 100644 --- a/include/ensmallen_bits/nsga3/nsga3_impl.hpp +++ b/include/ensmallen_bits/nsga3/nsga3_impl.hpp @@ -16,17 +16,17 @@ #include "nsga3.hpp" #include +#include "normalization.hpp" namespace ens { -template -inline NSGA3::NSGA3( - const arma::Cube& referencePoints, +template +inline NSGA3::NSGA3( + const arma::Mat& referencePoints, const size_t populationSize, const size_t maxGenerations, const double crossoverProb, const double distributionIndex, - const double epsilon, const double eta, const arma::vec& lowerBound, const arma::vec& upperBound): @@ -37,20 +37,18 @@ inline NSGA3::NSGA3( maxGenerations(maxGenerations), crossoverProb(crossoverProb), distributionIndex(distributionIndex), - epsilon(epsilon), eta(eta), lowerBound(lowerBound), upperBound(upperBound) { /* Nothing to do here. */ } -template -inline NSGA3::NSGA3( - const arma::Cube& referencePoints, +template +inline NSGA3::NSGA3( + const arma::Mat& referencePoints, const size_t populationSize, const size_t maxGenerations, const double crossoverProb, const double distributionIndex, - const double epsilon, const double eta, const double lowerBound, const double upperBound): @@ -61,17 +59,17 @@ inline NSGA3::NSGA3( maxGenerations(maxGenerations), crossoverProb(crossoverProb), distributionIndex(distributionIndex), - epsilon(epsilon), eta(eta), lowerBound(lowerBound * arma::ones(1, 1)), upperBound(upperBound * arma::ones(1, 1)) { /* Nothing to do here. */ } //! Optimize the function. -template -template +template -typename MatType::elem_type NSGA3::Optimize( +typename MatType::elem_type NSGA3::Optimize( std::tuple& objectives, MatType& iterateIn, CallbackTypes&&... callbacks) @@ -111,8 +109,11 @@ typename MatType::elem_type NSGA3::Optimize( numObjectives = sizeof...(ArbitraryFunctionType); numVariables = iterate.n_rows; + assert(numObjectives == referencePoints.n_rows && "The dimensions of " + "reference points do not match the number of functions."); + // Cache calculated objectives. - std::vector > calculatedObjectives(populationSize); + std::vector > calculatedObjectives(populationSize); // Population size reserved to 2 * populationSize + 1 to accommodate // for the size of intermediate candidate population. @@ -144,6 +145,7 @@ typename MatType::elem_type NSGA3::Optimize( // Constrain all genes to be within bounds. population[i] = arma::min(arma::max(population[i], castedLowerBound), castedUpperBound); } + Normalization hpn(numObjectives); Info << "NSGA3 initialized successfully. Optimization started." << std::endl; @@ -152,6 +154,7 @@ typename MatType::elem_type NSGA3::Optimize( for (size_t generation = 1; generation <= maxGenerations && !terminate; generation++) { + std::cout << generation << std::endl; // Create new population of candidate from the present elite population. // Have P_t, generate G_t using P_t. BinaryTournamentSelection(population, castedLowerBound, castedUpperBound); @@ -159,65 +162,89 @@ typename MatType::elem_type NSGA3::Optimize( // Evaluate the objectives for the new population. calculatedObjectives.resize(population.size()); std::fill(calculatedObjectives.begin(), calculatedObjectives.end(), - arma::Col(numObjectives, arma::fill::zeros)); + arma::Col(numObjectives, arma::fill::zeros)); + EvaluateObjectives(population, objectives, calculatedObjectives); // Perform fast non dominated sort on P_t ∪ G_t. ranks.resize(population.size()); - FastNonDominatedSort(fronts, ranks, calculatedObjectives); + FastNonDominatedSort>(fronts, ranks, calculatedObjectives); + + hpn.update(calculatedObjectives, fronts[0]); + arma::Col denom = hpn.NadirPoint() - hpn.IdealPoint(); // S_t and P_t+1 declared. - std::vector selectedPoints(fronts[0].begin(), fronts[0].end()); - std::vector nextPopulation(fronts[0].begin(), fronts[0].end()); + std::vector selectedPoints; + std::vector nextPopulation; size_t index = 0; - while (nextPopulation.size() + fronts.size() < populationSize) + while (nextPopulation.size() + fronts[index].size() < populationSize) { selectedPoints.insert(selectedPoints.end(), fronts[index].begin(), fronts[index].end()); nextPopulation.insert(nextPopulation.end(), fronts[index].begin(), fronts[index].end()); index++; } - selectedPoints.insert(selectedPoints.end(), fronts[++index].begin(), fronts[++index].end()); - size_t lastFront = index; - arma::Col idealPoint(calculatedObjectives[selectedPoints[0]]); - for (index = 1; index < selectedPoints.size(); index++) + if(nextPopulation.size() != populationSize) { - idealPoint = arma::min(idealPoint, calculatedObjectives[selectedPoints[index]]); - } + selectedPoints.insert(selectedPoints.end(), fronts[index].begin(), fronts[index].end()); - for (index = 0; index < selectedPoints.size(); index++) - { - calculatedObjectives[selectedPoints[index]] = - calculatedObjectives[selectedPoints[index]] - idealPoint; - } + size_t lastFront = index; - arma::Col normalize(numObjectives, arma::fill::zeros); - arma::Row extreme(numObjectives, arma::fill::zeros); - - // Find the extreme points and normalize S_t - FindExtremePoints>(extreme, calculatedObjectives, - selectedPoints); - NormalizeFront>(calculatedObjectives, normalize, - selectedPoints, extreme); - - // Find the associated reference directions to the selected points. - arma::urowvec refIndex(selectedPoints.size()); - arma::Row dists(selectedPoints.size()); - Associate>(refIndex, dists, calculatedObjectives, - selectedPoints); - - // Calculate the niche count of S_t and performing the niching operation. - arma::Row count(referencePoints.n_slices); - NicheCount(count, refIndex); - Niching(selectedPoints.size() - nextPopulation.size(), count, refIndex, - dists, fronts[lastFront], nextPopulation); - + for (index = 0; index < selectedPoints.size(); index++) + { + calculatedObjectives[selectedPoints[index]] = + calculatedObjectives[selectedPoints[index]] - hpn.IdealPoint(); + calculatedObjectives[selectedPoints[index]] = + calculatedObjectives[selectedPoints[index]] / denom; + } + + // Find the associated reference directions to the selected points. + arma::urowvec refIndex(selectedPoints.size()); + arma::Row dists(selectedPoints.size()); + /*std::cout << "["; + for (size_t i = 0; i < referencePoints.n_cols; i++) + { + std::cout << "["; + for (size_t j = 0; j < referencePoints.n_rows; j++) + { + std::cout << referencePoints(j, i) << ","; + } + std::cout << "],"; + } + std::cout << "]" << std::endl; + std::cout << "["; + for (size_t i = 0; i < selectedPoints.size(); i++) + { + std::cout << "["; + for (size_t j = 0; j < 3; j++) + { + std::cout << calculatedObjectives[selectedPoints[i]][j] << ","; + } + std::cout << "],"; + }*/ + //std::cout << "]" << std::endl; + Associate>(refIndex, dists, calculatedObjectives, + selectedPoints); + + /*std::cout << "["; + for (size_t i = 0; i < refIndex.n_elem; i++) + { + std::cout << refIndex[i] << ","; + }*/ + // Calculate the niche count of S_t and performing the niching operation. + arma::Row count(referencePoints.n_cols, arma::fill::zeros); + + NicheCount(count, refIndex, nextPopulation); + Niching(populationSize - nextPopulation.size(), count, refIndex, + dists, fronts[lastFront], nextPopulation); + } for (size_t i : nextPopulation) { tempPopulation.push_back(population[i]); } population = tempPopulation; + tempPopulation.erase(tempPopulation.begin(), tempPopulation.end()); terminate |= Callback::GenerationalStepTaken(*this, objectives, iterate, calculatedObjectives, fronts, callbacks...); @@ -253,7 +280,7 @@ typename MatType::elem_type NSGA3::Optimize( ElemType performance = std::numeric_limits::max(); - for (const arma::Col& objective: calculatedObjectives) + for (const arma::Col& objective: calculatedObjectives) if (arma::accu(objective) < performance) performance = arma::accu(objective); @@ -261,31 +288,31 @@ typename MatType::elem_type NSGA3::Optimize( } //! No objectives to evaluate. -template +template template typename std::enable_if::type -NSGA3::EvaluateObjectives( +NSGA3::EvaluateObjectives( std::vector&, std::tuple&, - std::vector&) + std::vector >&) { // Nothing to do here. } //! Evaluate the objectives for the entire population. -template +template template typename std::enable_if::type -NSGA3::EvaluateObjectives( +NSGA3::EvaluateObjectives( std::vector& population, std::tuple& objectives, - std::vector& calculatedObjectives) + std::vector >& calculatedObjectives) { - for (size_t i = 0; i < populationSize; i++) + for (size_t i = 0; i < population.size(); i++) { calculatedObjectives[i](I) = std::get(objectives).Evaluate(population[i]); EvaluateObjectives(population, objectives, @@ -294,15 +321,16 @@ NSGA3::EvaluateObjectives( } //! Reproduce and generate new candidates. +template template -inline void NSGA3::BinaryTournamentSelection( +inline void NSGA3::BinaryTournamentSelection( std::vector& population, const MatType& lowerBound, const MatType& upperBound) { std::vector children; - while (children.size() < population.size()) + while (children.size() < populationSize) { // Choose two random parents for reproduction from the elite population. size_t indexA = arma::randi(arma::distr_param(0, populationSize - 1)); @@ -339,13 +367,14 @@ inline void NSGA3::BinaryTournamentSelection( } //! Perform simulated binary crossover (SBX) of genes for the children. +template template -inline void NSGA3::Crossover(MatType& childA, - MatType& childB, - const MatType& parentA, - const MatType& parentB, - const MatType& lowerBound, - const MatType& upperBound) +inline void NSGA3::Crossover(MatType& childA, + MatType& childB, + const MatType& parentA, + const MatType& parentB, + const MatType& lowerBound, + const MatType& upperBound) { //! Generates a child from two parent individuals // according to the polynomial probability distribution. @@ -397,11 +426,12 @@ inline void NSGA3::Crossover(MatType& childA, } //! Perform Polynomial mutation of the candidate. +template template -inline void NSGA3::Mutate(MatType& candidate, - double mutationRate, - const MatType& lowerBound, - const MatType& upperBound) +inline void NSGA3::Mutate(MatType& candidate, + double mutationRate, + const MatType& lowerBound, + const MatType& upperBound) { const size_t numVariables = candidate.n_rows; for (size_t geneIdx = 0; geneIdx < numVariables; ++geneIdx) @@ -436,110 +466,35 @@ inline void NSGA3::Mutate(MatType& candidate, candidate = arma::min(arma::max(candidate, lowerBound), upperBound); } -template -template -inline void NSGA3::NormalizeFront( - std::vector& calculatedObjectives, - ColType& normalization, - const std::vector& front, - const arma::Row& extreme) -{ - arma::Mat vectorizedObjectives(numObjectives, front.size()); - for (size_t i = 0; i < front.size(); i++) - { - vectorizedObjectives.col(i) = calculatedObjectives[front[i]]; - } - - if(front.size() < numObjectives) - { - normalization = arma::max(vectorizedObjectives, 1); - return; - } - ColType temp; - arma::uvec unique = arma::find_unique(extreme); - if (extreme.n_elem != unique.n_elem) - { - normalization = arma::max(vectorizedObjectives, 1); - return; - } - ColType one(front.size(), arma::fill::ones); - ColType hyperplane = arma::solve( - vectorizedObjectives.t(), one); - if (hyperplane.has_inf() || hyperplane.has_nan() || (arma::accu(hyperplane < 0.0) > 0)) - { - normalization = arma::max(vectorizedObjectives, 1); - } - else - { - normalization = 1. / hyperplane; - if (hyperplane.has_inf() || hyperplane.has_nan()) - { - normalization = arma::max(vectorizedObjectives, 1); - } - } - normalization = normalization + (normalization == 0); -} - -//! Find the index of the of the extreme points in the given front. -template -template -void NSGA3::FindExtremePoints(arma::Row& indexes, - std::vector& calculatedObjectives, - const std::vector& front) -{ - typedef typename MatType::elem_type ElemType; - - if(numObjectives >= front.size()) - { - indexes = arma::linspace>(0, front.size() - 1, front.size()); - return; - } - - arma::Mat W(numObjectives, numObjectives, arma::fill::eye); - W = W + 1e-6; - std::vector selected(front.size()); - arma::Col z(numObjectives, arma::fill::zeros); - arma::Row dists; - for (size_t i = 0; i < numObjectives; i++) - { - PointToLineDistance(dists, calculatedObjectives, front, z, W.col(i)); - for (size_t j = 0; j < front.size(); j++) - if (selected[j]){dists[j] = arma::datum::inf;} - indexes[i] = dists.index_min(); - selected[dists.index_min()] = true; - } -} - //! Find the distance of a front from a line formed by two points. -template +template template -void NSGA3::PointToLineDistance( - arma::Row& distances, - std::vector& calculatedObjectives, +inline void NSGA3::PointToLineDistance( + arma::Row& distances, + const std::vector& calculatedObjectives, const std::vector& front, const ColType& pointA, const ColType& pointB) { - typedef typename MatType::elem_type ElemType; - arma::Row distancesTemp(front.size()); + arma::Row distancesTemp(front.size()); ColType ba = pointB - pointA; ColType pa; for (size_t i = 0; i < front.size(); i++) { - size_t ind = front[i]; + size_t ind = front[i]; - pa = (calculatedObjectives[ind] - pointA); - double t = arma::dot(pa, ba) / arma::dot(ba, ba); - distancesTemp[i] = arma::accu(arma::pow((pa - t * ba), 2)); + pa = (calculatedObjectives[ind] - pointA); + double t = arma::dot(pa, ba) / arma::dot(ba, ba); + distancesTemp[i] = std::pow(arma::accu(arma::pow((pa - t * ba), 2)), 0.5); } distances = distancesTemp; } //! Sort population into Pareto fronts. -template +template template -inline void NSGA3::FastNonDominatedSort( +inline void NSGA3::FastNonDominatedSort( std::vector >& fronts, std::vector& ranks, std::vector& calculatedObjectives) @@ -551,16 +506,16 @@ inline void NSGA3::FastNonDominatedSort( fronts.clear(); fronts.push_back(std::vector()); - for (size_t p = 0; p < populationSize; p++) + for (size_t p = 0; p < calculatedObjectives.size(); p++) { dominated[p] = std::set(); dominationCount[p] = 0; - for (size_t q = 0; q < populationSize; q++) + for (size_t q = 0; q < calculatedObjectives.size(); q++) { - if (Dominates(calculatedObjectives, p, q)) + if (Dominates>(calculatedObjectives, p, q)) dominated[p].insert(q); - else if (Dominates(calculatedObjectives, q, p)) + else if (Dominates>(calculatedObjectives, q, p)) dominationCount[p] += 1; } @@ -598,23 +553,24 @@ inline void NSGA3::FastNonDominatedSort( fronts.pop_back(); } -template -void NSGA3::Niching(size_t K, - arma::Row& nicheCount, - const arma::urowvec& refIndex, - const arma::Row& dists, - const std::vector& front, - std::vector& population) +template +inline void NSGA3::Niching(size_t K, + arma::Row& nicheCount, + const arma::urowvec& refIndex, + const arma::Row& dists, + const std::vector& front, + std::vector& population) { - arma::Row popMask(population.size()); + arma::Row popMask(front.size(), arma::fill::zeros); + int nextPopSize = population.size(); size_t k = 0; while (k < K) { size_t jMin = arma::index_min(nicheCount); std::vector I; - for (size_t i : front) + for (size_t i = 0; i < front.size(); i++) { - if(refIndex[i] == jMin && popMask[i]) + if(refIndex[nextPopSize + i] == jMin && !popMask[i]) I.push_back(i); } if (I.size() != 0) @@ -624,58 +580,62 @@ void NSGA3::Niching(size_t K, { for (size_t i = 0; i < I.size(); i++) { - if(dists[I[i]] < dists[I[min]]) + if(dists[nextPopSize + I[i]] < dists[nextPopSize + I[min]]) { min = i; } } } - population.push_back(I[min]); + population.push_back(front[I[min]]); nicheCount[jMin] += 1; - popMask[I[min]] = false; - k += 1; + popMask[I[min]] = 1; + k++; } else { - nicheCount[jMin] = arma::datum::inf; + nicheCount[jMin] = 100000; } } } -template +template template -void NSGA3::Associate(arma::urowvec& refIndex, - arma::Row& dists, - const std::vector& calculatedObjectives, - const std::vector& St) +inline void NSGA3::Associate( + arma::urowvec& refIndex, + arma::Row& dists, + const std::vector& calculatedObjectives, + const std::vector& St) { - arma::Mat d(refIndex.n_elem, St.size()); - ColType zero(arma::size(calculatedObjectives),arma::fill::zeros); - for (size_t i = 0; i < referencePoints.n_slices; i++) + arma::Mat d(referencePoints.n_cols, St.size()); + ColType zero(arma::size(calculatedObjectives[0]),arma::fill::zeros); + arma::Row temp; + for (size_t i = 0; i < referencePoints.n_cols; i++) { - PointToLineDistance(d.row(i), calculatedObjectives, St, - referencePoints.slice(i), zero); + PointToLineDistance(temp, calculatedObjectives, St, + zero, referencePoints.col(i)); + d.row(i) = temp; } refIndex = arma::index_min(d, 0); dists = arma::min(d, 0); } -template -void NSGA3::NicheCount(arma::Row& count, - const arma::urowvec& refIndex) +template +inline void NSGA3::NicheCount(arma::Row& count, + const arma::urowvec& refIndex, + const std::vector& nextPopulation) { - for (size_t i = 0; i < refIndex.n_elem; i++) + for (size_t i = 0; i < nextPopulation.size(); i++) { count[refIndex[i]] += 1; } } //! Check if a candidate Pareto dominates another candidate. -template +template template -inline bool NSGA3::Dominates( +inline bool NSGA3::Dominates( std::vector& calculatedObjectives, size_t candidateP, size_t candidateQ) diff --git a/include/ensmallen_bits/problems/dtlz/dtlz1_function.hpp b/include/ensmallen_bits/problems/dtlz/dtlz1_function.hpp index d1ff5223c..46b815de4 100644 --- a/include/ensmallen_bits/problems/dtlz/dtlz1_function.hpp +++ b/include/ensmallen_bits/problems/dtlz/dtlz1_function.hpp @@ -179,10 +179,6 @@ namespace test { { value = value * (1. - coords[stop]); } - else - { - value = value * coords[stop]; - } value = value * (1. + dtlz.g(coords)[0]); return value; From 9e8208adfcf4cbf6a8f38b9caaca6f2d98957c55 Mon Sep 17 00:00:00 2001 From: IWNMWE Date: Thu, 22 Aug 2024 02:09:59 +0530 Subject: [PATCH 06/11] remove cout --- include/ensmallen_bits/nsga3/nsga3_impl.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/ensmallen_bits/nsga3/nsga3_impl.hpp b/include/ensmallen_bits/nsga3/nsga3_impl.hpp index ebe513509..6e4886c80 100644 --- a/include/ensmallen_bits/nsga3/nsga3_impl.hpp +++ b/include/ensmallen_bits/nsga3/nsga3_impl.hpp @@ -154,7 +154,6 @@ typename MatType::elem_type NSGA3::Optimize( for (size_t generation = 1; generation <= maxGenerations && !terminate; generation++) { - std::cout << generation << std::endl; // Create new population of candidate from the present elite population. // Have P_t, generate G_t using P_t. BinaryTournamentSelection(population, castedLowerBound, castedUpperBound); From 263d160dec659b9adab098805578dd977f6bde93 Mon Sep 17 00:00:00 2001 From: IWNMWE Date: Thu, 22 Aug 2024 02:22:38 +0530 Subject: [PATCH 07/11] Remove cout --- include/ensmallen_bits/nsga3/nsga3_impl.hpp | 28 +-------------------- 1 file changed, 1 insertion(+), 27 deletions(-) diff --git a/include/ensmallen_bits/nsga3/nsga3_impl.hpp b/include/ensmallen_bits/nsga3/nsga3_impl.hpp index 6e4886c80..a3d3b19f9 100644 --- a/include/ensmallen_bits/nsga3/nsga3_impl.hpp +++ b/include/ensmallen_bits/nsga3/nsga3_impl.hpp @@ -201,36 +201,10 @@ typename MatType::elem_type NSGA3::Optimize( // Find the associated reference directions to the selected points. arma::urowvec refIndex(selectedPoints.size()); arma::Row dists(selectedPoints.size()); - /*std::cout << "["; - for (size_t i = 0; i < referencePoints.n_cols; i++) - { - std::cout << "["; - for (size_t j = 0; j < referencePoints.n_rows; j++) - { - std::cout << referencePoints(j, i) << ","; - } - std::cout << "],"; - } - std::cout << "]" << std::endl; - std::cout << "["; - for (size_t i = 0; i < selectedPoints.size(); i++) - { - std::cout << "["; - for (size_t j = 0; j < 3; j++) - { - std::cout << calculatedObjectives[selectedPoints[i]][j] << ","; - } - std::cout << "],"; - }*/ - //std::cout << "]" << std::endl; + Associate>(refIndex, dists, calculatedObjectives, selectedPoints); - /*std::cout << "["; - for (size_t i = 0; i < refIndex.n_elem; i++) - { - std::cout << refIndex[i] << ","; - }*/ // Calculate the niche count of S_t and performing the niching operation. arma::Row count(referencePoints.n_cols, arma::fill::zeros); From 3c3a8dd8a00ce9cbef75ab4f807d999fa2bf889b Mon Sep 17 00:00:00 2001 From: IWNMWE Date: Fri, 27 Sep 2024 00:14:15 +0530 Subject: [PATCH 08/11] change evaluate --- include/ensmallen_bits/nsga3/nsga3.hpp | 1 + include/ensmallen_bits/nsga3/nsga3_impl.hpp | 20 +++++++++++--------- 2 files changed, 12 insertions(+), 9 deletions(-) diff --git a/include/ensmallen_bits/nsga3/nsga3.hpp b/include/ensmallen_bits/nsga3/nsga3.hpp index d37fe1715..9a58f025f 100644 --- a/include/ensmallen_bits/nsga3/nsga3.hpp +++ b/include/ensmallen_bits/nsga3/nsga3.hpp @@ -202,6 +202,7 @@ class NSGA3 return rcFront; } + private: /** * Evaluate objectives for the elite population. * diff --git a/include/ensmallen_bits/nsga3/nsga3_impl.hpp b/include/ensmallen_bits/nsga3/nsga3_impl.hpp index a3d3b19f9..c20a3cc74 100644 --- a/include/ensmallen_bits/nsga3/nsga3_impl.hpp +++ b/include/ensmallen_bits/nsga3/nsga3_impl.hpp @@ -156,6 +156,7 @@ typename MatType::elem_type NSGA3::Optimize( { // Create new population of candidate from the present elite population. // Have P_t, generate G_t using P_t. + std::cout << generation << std::endl; BinaryTournamentSelection(population, castedLowerBound, castedUpperBound); // Evaluate the objectives for the new population. @@ -224,22 +225,23 @@ typename MatType::elem_type NSGA3::Optimize( } EvaluateObjectives(population, objectives, calculatedObjectives); // Set the candidates from the Pareto Set as the output. - paretoSet.set_size(population[0].n_rows, population[0].n_cols, fronts[0].size()); + paretoSet.set_size(population[0].n_rows, population[0].n_cols, + population.size()); // The Pareto Set is stored, can be obtained via ParetoSet() getter. - for (size_t solutionIdx = 0; solutionIdx < fronts[0].size(); ++solutionIdx) + for (size_t solutionIdx = 0; solutionIdx < population.size(); ++solutionIdx) { paretoSet.slice(solutionIdx) = - arma::conv_to::from(population[fronts[0][solutionIdx]]); + arma::conv_to::from(population[solutionIdx]); } // Set the candidates from the Pareto Front as the output. - paretoFront.set_size(calculatedObjectives[0].n_rows, calculatedObjectives[0].n_cols, - fronts[0].size()); + paretoFront.set_size(calculatedObjectives[0].n_rows, + calculatedObjectives[0].n_cols, population.size()); // The Pareto Front is stored, can be obtained via ParetoFront() getter. - for (size_t solutionIdx = 0; solutionIdx < fronts[0].size(); ++solutionIdx) + for (size_t solutionIdx = 0; solutionIdx < population.size(); ++solutionIdx) { paretoFront.slice(solutionIdx) = - arma::conv_to::from(calculatedObjectives[fronts[0][solutionIdx]]); + arma::conv_to::from(calculatedObjectives[solutionIdx]); } // Clear rcFront, in case it is later requested by the user for reverse @@ -288,9 +290,9 @@ NSGA3::EvaluateObjectives( for (size_t i = 0; i < population.size(); i++) { calculatedObjectives[i](I) = std::get(objectives).Evaluate(population[i]); - EvaluateObjectives(population, objectives, - calculatedObjectives); } + EvaluateObjectives(population, objectives, + calculatedObjectives); } //! Reproduce and generate new candidates. From 40142e5644bbb11b101f1cecb483f5e4e0196f9a Mon Sep 17 00:00:00 2001 From: IWNMWE Date: Fri, 27 Sep 2024 00:15:43 +0530 Subject: [PATCH 09/11] remove cout --- include/ensmallen_bits/nsga3/nsga3_impl.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/include/ensmallen_bits/nsga3/nsga3_impl.hpp b/include/ensmallen_bits/nsga3/nsga3_impl.hpp index c20a3cc74..9b0e5c1d9 100644 --- a/include/ensmallen_bits/nsga3/nsga3_impl.hpp +++ b/include/ensmallen_bits/nsga3/nsga3_impl.hpp @@ -156,7 +156,6 @@ typename MatType::elem_type NSGA3::Optimize( { // Create new population of candidate from the present elite population. // Have P_t, generate G_t using P_t. - std::cout << generation << std::endl; BinaryTournamentSelection(population, castedLowerBound, castedUpperBound); // Evaluate the objectives for the new population. From 7a2ec9269330320e9b8ff6869ea1d5aaaed7e0f1 Mon Sep 17 00:00:00 2001 From: IWNMWE Date: Sat, 30 Nov 2024 10:47:43 +0530 Subject: [PATCH 10/11] nsga3: typo fix --- include/ensmallen_bits/nsga3/nsga3.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/ensmallen_bits/nsga3/nsga3.hpp b/include/ensmallen_bits/nsga3/nsga3.hpp index 9a58f025f..70d843441 100644 --- a/include/ensmallen_bits/nsga3/nsga3.hpp +++ b/include/ensmallen_bits/nsga3/nsga3.hpp @@ -252,7 +252,7 @@ class NSGA3 * @param parentA First parent from elite population. * @param parentB Second parent from elite population. * @param lowerBound Lower Bound of the offspring. - * @param upperBound Upper Boundn of the offspring. + * @param upperBound Upper Bound of the offspring. */ template void Crossover(MatType& childA, From db9c81b925e1e50c8b0e1f34bd66a792f16204f2 Mon Sep 17 00:00:00 2001 From: IWNMWE Date: Sat, 30 Nov 2024 10:49:08 +0530 Subject: [PATCH 11/11] nsga3: documentation added --- doc/optimizers.md | 69 ++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 68 insertions(+), 1 deletion(-) diff --git a/doc/optimizers.md b/doc/optimizers.md index b5f31c470..c739cb2c1 100644 --- a/doc/optimizers.md +++ b/doc/optimizers.md @@ -2214,7 +2214,7 @@ size equal to that of the starting population. | `double`, `arma::vec` | **`lowerBound`** | Lower bound of the coordinates on the coordinates of the whole population during the search process. | `0` | | `double`, `arma::vec` | **`upperBound`** | Lower bound of the coordinates on the coordinates of the whole population during the search process. | `1` | -Note that the parameters `lowerBound` and `upperBound` are overloaded. Data types of `double` or `arma::mat` may be used. If they are initialized as single values of `double`, then the same value of the bound applies to all the axes, resulting in an initialization following a uniform distribution in a hypercube. If they are initialized as matrices of `arma::mat`, then the value of `lowerBound[i]` applies to axis `[i]`; similarly, for values in `upperBound`. This results in an initialization following a uniform distribution in a hyperrectangle within the specified bounds. +Note that the parameters `lowerBound` and `upperBound` are overloaded. Data types of `double` or `arma::mat` may be used. If they are initialized as single values of `double`, then the same value of the bound applies to all the axes, resulting in an initialization following a uniform distribution in a hypercube. If they are initialized as matrices of `arma::mat`, then the value of `lowerBound[i]` applies to axis `[i]`; similarly, for values in `upperBound`. This results in an initialization following a uniform distribution in a hyperrectangle within the specified bounds.The `ElementType` parameter refer Attributes of the optimizer may also be changed via the member methods `PopulationSize()`, `MaxGenerations()`, `CrossoverRate()`, `MutationProbability()`, `MutationStrength()`, `Epsilon()`, `LowerBound()` and `UpperBound()`. @@ -2252,6 +2252,73 @@ arma::cube bestFront = opt.ParetoFront(); * [Multi-objective functions](#multi-objective-functions) * [Performance Indicators](#performance-indicators) +## NSGA3 + +*An optimizer for arbitrary multi-objective functions.*\ + +NSGA-III (Non-dominated Sorting Genetic Algorithm III) is a multi-objective optimization algorithm designed to handle problems with many objectives more effectively than its predecessor, NSGA-II. Like NSGA-II, it starts by generating an initial population and creates a new population of offspring at each iteration. The combined population (parents and offspring) is sorted into non-dominated fronts. However, NSGA-III introduces a reference-point-based selection mechanism, ensuring a well-distributed set of solutions across the Pareto front. This approach improves performance for higher-dimensional objective spaces by maintaining diversity and focusing on convergence. + +#### Constructors + + * `NSGA3()` + * `NSGA3(`_`referencePoints, populationSize, maxGenerations, crossoverProb, distributionIndex, eta, lowerBound, upperBound`_`)` + + The _`ElementType`_ template parameter refers to the type of elements used in + initiallization of matrices. + +#### Attributes + +| **type** | **name** | **description** | **default** | +|----------|----------|-----------------|-------------| +| `size_t` | **`populationSize`** | The number of candidates in the population. This should be at least 4 in size and a multiple of 4. | `100` | +| `size_t` | **`maxGenerations`** | The maximum number of generations allowed for NSGA3. | `2000` | +| `double` | **`crossoverProb`** | Probability that a crossover will occur. | `0.6` | +| `double` | **`distributionIndex`** | The crowding degree of the mutation. | `20` | +| `double` | **`eta`** | The distance parameters of the crossover distribution. | `20` | +| `double`, `arma::vec` | **`lowerBound`** | Lower bound of the coordinates on the coordinates of the whole population during the search process. | `0` | +| `double`, `arma::vec` | **`upperBound`** | Lower bound of the coordinates on the coordinates of the whole population during the search process. | `1` | + +Note that the parameters `lowerBound` and `upperBound` are overloaded. Data types of `double` or `arma::mat` may be used. If they are initialized as single values of `double`, then the same value of the bound applies to all the axes, resulting in an initialization following a uniform distribution in a hypercube. If they are initialized as matrices of `arma::mat`, then the value of `lowerBound[i]` applies to axis `[i]`; similarly, for values in `upperBound`. This results in an initialization following a uniform distribution in a hyperrectangle within the specified bounds. + +Attributes of the optimizer may also be changed via the member methods +`PopulationSize()`, `MaxGenerations()`, `CrossoverRate()`, `DistributionIndex()`, `Eta()`, `ReferencePoints()`, `LowerBound()` and `UpperBound()`. + +#### Examples: + +
+Click to collapse/expand example code. + + +```c++ +DTLZ1 DtlzOne; +const double lowerBound = 0; +const double upperBound = 1; +const double expectedLowerBound = 0.5; +const double expectedUpperBound = 0.5; + +Uniform generator; +arma::Cube refPoints(3,1,136); +arma::mat g = generator.Generate(3,91,0); + +NSGA3 opt(g, 120, 100, 1, 20, 30, 0, 1); + +typedef decltype(DtlzOne.objectiveF1) ObjectiveTypeA; +typedef decltype(DtlzOne.objectiveF2) ObjectiveTypeB; +typedef decltype(DtlzOne.objectiveF3) ObjectiveTypeC; + +arma::Col coords = DtlzOne.GetInitialPoint(); +std::tuple objectives = + DtlzOne.GetObjectives(); +opt.Optimize(objectives, coords); +``` + +
+ +#### See also: + + * [NSGA-III Algorithm](https://www.egr.msu.edu/~kdeb/papers/k2012009.pdf) + * [Reference point based multi-objective optimization using evolutionary algorithms](https://www.egr.msu.edu/~kdeb/papers/k2005012.pdf) + ## OptimisticAdam *An optimizer for [differentiable separable functions](#differentiable-separable-functions).*