diff --git a/doc/function_types.md b/doc/function_types.md index 418ac3501..b5f282136 100644 --- a/doc/function_types.md +++ b/doc/function_types.md @@ -884,6 +884,86 @@ arma::cube bestFront = optimizer.ParetoFront(); +### Performance Indicators + +Performance indicators in multiobjective optimization provide essential metrics +for evaluating solution quality, such as convergence to the Pareto front and solution +diversity. + +The ensmallen library offers three such indicators, aiding in the assessment and comparison +of different optimization methods: + +#### Epsilon + +Epsilon metric is a performance metric used in multi-objective optimization which measures +the smallest factor by which a set of solution objectives must be scaled to dominate a +reference set of solutions. Specifically, given a set of Pareto-optimal solutions, the +epsilon indicator finds the minimum value ϵ such that each solution in the set is at +least as good as every solution in the reference set when the objectives are scaled by ϵ. + +
+Click to collapse/expand example code. + + +```c++ +arma::cube referenceFront(2, 1, 3); +double tol = 1e-10; +referenceFront.slice(0) = arma::vec{0.01010101, 0.89949622}; +referenceFront.slice(1) = arma::vec{0.02020202, 0.85786619}; +referenceFront.slice(2) = arma::vec{0.03030303, 0.82592234}; +arma::cube front = referenceFront * 1.1; +// eps is approximately 1.1 +double eps = Epsilon::Evaluate(front, referenceFront); +``` +
+ +#### IGD + +Inverse Generational Distance (IGD) is a performance metric used in multi-objective optimization +to evaluate the quality of a set of solutions relative to a reference set, typically representing +the true Pareto front. IGD measures the average distance from each point in the reference set to +the closest point in the obtained solution set. + +
+Click to collapse/expand example code. + + +```c++ +arma::cube referenceFront(2, 1, 3); +double tol = 1e-10; +referenceFront.slice(0) = arma::vec{0.01010101, 0.89949622}; +referenceFront.slice(1) = arma::vec{0.02020202, 0.85786619}; +referenceFront.slice(2) = arma::vec{0.03030303, 0.82592234}; +arma::cube front = referenceFront * 1.1; +// The third parameter is the power constant in the distance formula. +// IGD is approximately 0.05329 +double igd = IGD::Evaluate(front, referenceFront, 1); +``` +
+ +#### IGD Plus + +IGD Plus (IGD+) is a variant of the Inverse Generational Distance (IGD) metric used in multi-objective +optimization. It refines the traditional IGD metric by incorporating a preference for Pareto-dominance +in the distance calculation. This modification helps IGD+ better reflect both the convergence to the +Pareto front and the diversity of the solution set. + +
+Click to collapse/expand example code. + + +```c++ +arma::cube referenceFront(2, 1, 3); +double tol = 1e-10; +referenceFront.slice(0) = arma::vec{0.01010101, 0.89949622}; +referenceFront.slice(1) = arma::vec{0.02020202, 0.85786619}; +referenceFront.slice(2) = arma::vec{0.03030303, 0.82592234}; +arma::cube front = referenceFront * 1.1; +// IGDPlus is approximately 0.05329 +double igdPlus = IGDPlus::Evaluate(front, referenceFront); +``` +
+ *Note*: all multi-objective function optimizers have both the function `Optimize()` to find the best front, and also the function `ParetoFront()` to return all sets of solutions that are on the front. @@ -891,6 +971,11 @@ front. The following optimizers can be used with multi-objective functions: - [NSGA2](#nsga2) - [MOEA/D-DE](#moead) +- [AGEMOEA](#agemoea) + +#### See also: +* [Performance Assessment of Multiobjective Optimizers: An Analysis and Review](https://sop.tik.ee.ethz.ch/publicationListFiles/ztlf2003a.pdf) +* [Modified Distance Calculation in Generational Distance and Inverted Generational Distance](https://link.springer.com/chapter/10.1007/978-3-319-15892-1_8) ## Constrained functions diff --git a/doc/optimizers.md b/doc/optimizers.md index bb8a25073..b5f31c470 100644 --- a/doc/optimizers.md +++ b/doc/optimizers.md @@ -501,6 +501,97 @@ optimizer.Optimize(f, coordinates); * [Adam: A Method for Stochastic Optimization](http://arxiv.org/abs/1412.6980) (see section 7) * [Differentiable separable functions](#differentiable-separable-functions) +## AGEMOEA + +*An optimizer for arbitrary multi-objective functions.* + +Adaptive Geometry Estimation based Multi-Objective Evolutionary Algorithm (AGE-MOEA) is an optimization framework based on NSGA-II yet differs from it in replacing the crowding distance of NSGA-II by a survival score, for which calculations need the diversity and proximity of non-dominated sets. To simplify the computation of the survival score, in each generation, the geometry of the initial non-dominated subset is estimated by AGE-MOEA afterwards, this estimation which gets more accurate as the algorithm matures, is used as the geometry of the Pareto set. + +#### Constructors + + * `AGEMOEA()` + * `AGEMOEA(`_`populationSize, maxGenerations, crossoverProb, distributionIndex, epsilon, eta, lowerBound, upperBound`_`)` + +#### 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 AGEMOEA. | `2000` | +| `double` | **`crossoverProb`** | Probability that a crossover will occur. | `0.6` | +| `double` | **`distributionIndex`** | The crowding degree of the mutation. | `20` | +| `double` | **`epsilon`** | The value used internally to evaluate approximate equality in crowding distance based sorting. | `1e-6` | +| `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()`, `Epsilon()`, `LowerBound()` and `UpperBound()`. + +#### Examples + +
+Click to collapse/expand example code. + + +```c++ +SchafferFunctionN1 SCH; +arma::vec lowerBound("-1000"); +arma::vec upperBound("1000"); +AGEMOEA opt(50, 1000, 0.6, 20, 1e-6, 20, lowerBound, upperBound); + +typedef decltype(SCH.objectiveA) ObjectiveTypeA; +typedef decltype(SCH.objectiveB) ObjectiveTypeB; + +arma::mat coords = SCH.GetInitialPoint(); +std::tuple objectives = SCH.GetObjectives(); + +// obj will contain the minimum sum of objectiveA and objectiveB found on the best front. +double obj = opt.Optimize(objectives, coords); +// Now obtain the best front. +arma::cube bestFront = opt.ParetoFront(); +``` + +
+ +
+Click to collapse/expand example code. + + +```c++ +ZDT3<> ZDT_THREE(300); +const double lowerBound = 0; +const double upperBound = 1; + +AGEMOEA opt(50, 500, 0.8, 20, 1e-6, 20, lowerBound, upperBound); +typedef decltype(ZDT_THREE.objectiveF1) ObjectiveTypeA; +typedef decltype(ZDT_THREE.objectiveF2) ObjectiveTypeB; +bool success = true; +arma::mat coords = ZDT_THREE.GetInitialPoint(); +std::tuple objectives = ZDT_THREE.GetObjectives(); +opt.Optimize(objectives, coords); +const arma::cube bestFront = opt.ParetoFront(); + +NSGA2 opt2(50, 5000, 0.5, 0.5, 1e-3, 1e-6, lowerBound, upperBound); +// obj2 will contain the minimum sum of objectiveA and objectiveB found on the best front. +double obj2 = opt2.Optimize(objectives, coords); + +arma::cube NSGAFront = opt2.ParetoFront(); +// Get the IGD score for NSGA front using AGEMOEA as reference. +double igd = IGD::Evaluate(NSGAFront, bestFront, 1); +std::cout << igd << std::endl; +``` + +
+ +#### See also: + + * [An adaptive evolutionary algorithm based on non-euclidean geometry for many-objective optimization](https://doi.org/10.1145/3321707.3321839) + * [Multi-Objective Optimization in Wikipedia](https://en.wikipedia.org/wiki/Multi-objective_optimization) + * [Performance Indicators](#performance-indicators) + ## AMSBound *An optimizer for [differentiable separable functions](#differentiable-separable-functions).* @@ -2091,6 +2182,7 @@ arma::cube bestFront = opt.ParetoFront(); * [MOEA/D-DE Algorithm](https://ieeexplore.ieee.org/document/4633340) * [Multi-objective Functions in Wikipedia](https://en.wikipedia.org/wiki/Test_functions_for_optimization#Test_functions_for_multi-objective_optimization) * [Multi-objective functions](#multi-objective-functions) +* [Performance Indicators](#performance-indicators) ## NSGA2 @@ -2157,7 +2249,8 @@ arma::cube bestFront = opt.ParetoFront(); * [NSGA-II Algorithm](https://www.iitk.ac.in/kangal/Deb_NSGA-II.pdf) * [Multi-objective Functions in Wikipedia](https://en.wikipedia.org/wiki/Test_functions_for_optimization#Test_functions_for_multi-objective_optimization) - * [Multi-objective functions](#multi-objective-functions) + * [Multi-objective functions](#multi-objective-functions) + * [Performance Indicators](#performance-indicators) ## OptimisticAdam diff --git a/include/ensmallen.hpp b/include/ensmallen.hpp index bacc21026..c1e7f3984 100644 --- a/include/ensmallen.hpp +++ b/include/ensmallen.hpp @@ -66,6 +66,7 @@ #include "ensmallen_bits/utility/any.hpp" #include "ensmallen_bits/utility/arma_traits.hpp" #include "ensmallen_bits/utility/indicators/epsilon.hpp" +#include "ensmallen_bits/utility/indicators/igd.hpp" #include "ensmallen_bits/utility/indicators/igd_plus.hpp" // Contains traits, must be placed before report callback. @@ -111,6 +112,7 @@ #include "ensmallen_bits/katyusha/katyusha.hpp" #include "ensmallen_bits/lbfgs/lbfgs.hpp" #include "ensmallen_bits/lookahead/lookahead.hpp" +#include "ensmallen_bits/agemoea/agemoea.hpp" #include "ensmallen_bits/moead/moead.hpp" #include "ensmallen_bits/nsga2/nsga2.hpp" #include "ensmallen_bits/padam/padam.hpp" diff --git a/include/ensmallen_bits/agemoea/agemoea.hpp b/include/ensmallen_bits/agemoea/agemoea.hpp new file mode 100644 index 000000000..f7c5c2f79 --- /dev/null +++ b/include/ensmallen_bits/agemoea/agemoea.hpp @@ -0,0 +1,496 @@ +/** + * @file agemoea.hpp + * @author Satyam Shukla + * + * AGE-MOEA is a multi-objective optimization algorithm, widely used in + * many real-world applications. AGE-MOEA generates offsprings using + * crossover and mutation and then selects the next generation according + * to non-dominated-sorting and survival score comparison. + * + * 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_AGEMOEA_AGEMOEA_HPP +#define ENSMALLEN_AGEMOEA_AGEMOEA_HPP + +namespace ens { + +/** + * This class implements the AGEMOEA 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 + * @inproceedings{panichella2019adaptive, + * title={An adaptive evolutionary algorithm based on non-euclidean geometry for many-objective optimization}, + * author={Panichella, Annibale}, + * booktitle={Proceedings of the genetic and evolutionary computation conference}, + * pages={595--603}, + * year={2019} + * } + * @endcode + * + */ +class AGEMOEA +{ + public: + /** + * Constructor for the AGE-MOEA 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 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 distributionIndex The crowding degree of the mutation. + * @param epsilon The minimum difference required to distinguish between + * candidate solutions. + * @param eta The distance parameters of the crossover distribution. + * @param lowerBound Lower bound of the coordinates of the initial population. + * @param upperBound Upper bound of the coordinates of the initial population. + */ + AGEMOEA(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)); + + /** + * Constructor for the AGE-MOEA 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 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 distributionIndex The crowding degree of the mutation. + * @param epsilon The minimum difference required to distinguish between + * candidate solutions. + * @param eta The distance parameters of the crossover distribution + * @param lowerBound Lower bound of the coordinates of the initial population. + * @param upperBound Upper bound of the coordinates of the initial population. + */ + AGEMOEA(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); + + /** + * 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 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. + * @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; } + + //! Retrieve value of the distribution index. + double DistributionIndex() const { return distributionIndex; } + //! Modify the value of the distribution index. + double& DistributionIndex() { return distributionIndex; } + + //! 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; } + + //! 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. + */ + 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 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. + * @param lowerBound The lower bound of the objectives. + * @param upperBound The upper bound of the objectives. + */ + template + void Crossover(MatType& childA, + MatType& childB, + const MatType& parentA, + const MatType& parentB, + const MatType& lowerBound, + const MatType& upperBound); + + /** + * Mutate the coordinates for a candidate. + * + * @tparam MatType Type of matrix to optimize. + * @param candidate The candidate whose coordinates are being modified. + * @param mutationRate The probablity of a mutation to occur. + * @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, + 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); + + /** + * 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); + + /** + * Assigns Survival Score metric for sorting. + * + * @param front The previously generated Pareto fronts. + * @param idealPoint The ideal point of teh first front. + * @param calculatedObjectives The previously calculated objectives. + * @param survivalScore The Survival Score vector to be updated for each individual in the population. + * @param normalize The normlization vector of the fronts. + * @param dimension The dimension of the first front. + * @param fNum teh current front index. + */ + template + void SurvivalScoreAssignment( + const std::vector& front, + const arma::Col& idealPoint, + std::vector>& calculatedObjectives, + std::vector& survivalScore, + arma::Col& normalize, + double& dimension, + size_t fNum); + + /** + * The operator used in the AGE-MOEA survival score based sorting. + * + * If a candidate has a lower rank then it is preferred. + * Otherwise, if the ranks are equal then the candidate with the larger + * Survival Score 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 survivalScore The Survival score for each individual in + * the population. + * @return true if the first candidate is preferred, otherwise, false. + */ + template + bool SurvivalScoreOperator( + size_t idxP, + size_t idxQ, + const std::vector& ranks, + const std::vector& survivalScore); + + /** + * 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, + arma::Col& normalization, + const std::vector& front, + const arma::Row& extreme); + + /** + * Get the geometry information p of Lp norm (p > 0). + * + * @param calculatedObjectives The current population evaluated objectives. + * @param front The previously generated Pareto fronts. + * @param extreme The indexes of the extreme points in the front. + * @return The variable p in the Lp norm that best fits the geometry of the current front. + */ + template + double GetGeometry( + std::vector >& calculatedObjectives, + const std::vector& front, + const arma::Row& extreme); + + /** + * Finds the pairwise Lp distance between all the points in the front. + * + * @param final The current population evaluated objectives. + * @param calculatedObjectives The current population evaluated objectives. + * @param front The front of the current generation. + * @param dimension The calculated dimension of the front. + */ + template + void PairwiseDistance( + MatType& final, + std::vector >& calculatedObjectives, + const std::vector& front, + double dimension); + + /** + * 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 + * 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. + */ + template + void PointToLineDistance( + arma::Row& distances, + std::vector >& calculatedObjectives, + const std::vector& front, + const arma::Col& pointA, + const arma::Col& pointB); + + /** + * Find the Diversity score corresponding the solution S using the selected set. + * + * @param selected The current selected set. + * @param pairwiseDistance The current pairwise distance for the whole front. + * @param S The relative index of S being considered within the front. + * @return The diversity score for S which the sum of the two smallest elements. + */ + template + typename MatType::elem_type DiversityScore(std::set& selected, + const MatType& pairwiseDistance, + size_t S); + + //! 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 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; + + //! 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; + + //! 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 "agemoea_impl.hpp" + +#endif diff --git a/include/ensmallen_bits/agemoea/agemoea_impl.hpp b/include/ensmallen_bits/agemoea/agemoea_impl.hpp new file mode 100644 index 000000000..c643a7185 --- /dev/null +++ b/include/ensmallen_bits/agemoea/agemoea_impl.hpp @@ -0,0 +1,807 @@ +/** + * @file agemoea_impl.hpp + * @author Satyam Shukla + * + * Implementation of the AGEMOEA 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_AGEMOEA_AGEMOEA_IMPL_HPP +#define ENSMALLEN_AGEMOEA_AGEMOEA_IMPL_HPP + +#include "agemoea.hpp" +#include + +namespace ens { + +inline AGEMOEA::AGEMOEA(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) : + numObjectives(0), + numVariables(0), + populationSize(populationSize), + maxGenerations(maxGenerations), + crossoverProb(crossoverProb), + distributionIndex(distributionIndex), + epsilon(epsilon), + eta(eta), + lowerBound(lowerBound), + upperBound(upperBound) +{ /* Nothing to do here. */ } + +inline AGEMOEA::AGEMOEA(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) : + 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 +typename MatType::elem_type AGEMOEA::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("AGEMOEA::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 << "AGEMOEA 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); + + arma::Col idealPoint(calculatedObjectives[fronts[0][0]]); + for (size_t index = 1; index < fronts[0].size(); index++) + { + idealPoint = arma::min(idealPoint, + calculatedObjectives[fronts[0][index]]); + } + + // Perform survival score assignment. + survivalScore.resize(population.size()); + std::fill(survivalScore.begin(), survivalScore.end(), 0.); + double dimension; + arma::Col normalize(numObjectives, + arma::fill::zeros); + for (size_t fNum = 0; fNum < fronts.size(); fNum++) + { + SurvivalScoreAssignment(fronts[fNum], idealPoint, + calculatedObjectives, survivalScore, normalize, dimension, fNum); + } + + // 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); +} + +template +inline double AGEMOEA::GetGeometry( + std::vector >& calculatedObjectives, + const std::vector& front, + const arma::Row& extreme) +{ + arma::Row d; + arma::Col zero(numObjectives, arma::fill::zeros); + arma::Col one(numObjectives, arma::fill::ones); + + PointToLineDistance (d, calculatedObjectives, front, zero, one); + + for (size_t i = 0; i < extreme.size(); i++) + { + d[extreme[i]] = arma::datum::inf; + } + size_t index = arma::index_min(d); + double avg = arma::accu(calculatedObjectives[front[index]]) / static_cast (numObjectives); + double p = std::log(numObjectives) / std::log(1.0 / avg); + if (p <= 0.1 || std::isnan(p)) + p = 1.0; + + return p; +} + +//! Pairwise distance for each point in the given front. +template +inline void AGEMOEA::PairwiseDistance( + MatType& final, + std::vector >& calculatedObjectives, + const std::vector& front, + double dimension) +{ + for (size_t i = 0; i < front.size(); i++) + { + for (size_t j = i + 1; j < front.size(); j++) + { + final(i, j) = std::pow(arma::accu(arma::pow(arma::abs(calculatedObjectives[front[i]] - calculatedObjectives[front[j]]), dimension)), 1.0 / dimension); + final(j, i) = final(i, j); + } + } +} + +//! Find the index of the of the extreme points in the given front. +template +void AGEMOEA::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 AGEMOEA::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; +} + +//! Assign diversity score for a given point and teh selected set. +template +inline typename MatType::elem_type AGEMOEA::DiversityScore( + std::set& selected, + const MatType& pairwiseDistance, + size_t S) +{ + typedef typename MatType::elem_type ElemType; + ElemType m = arma::datum::inf; + ElemType m1 = arma::datum::inf; + std::set::iterator it; + for (it = selected.begin(); it != selected.end(); it++) + { + if (*it == S){ continue; } + if (pairwiseDistance(S, *it) < m) + { + m1 = m; + m = pairwiseDistance(S, *it); + } + else if (pairwiseDistance(S, *it) < m1) + { + m1 = pairwiseDistance(S, *it); + } + } + return m + m1; +} + +//! Assign survival score for a front of the population. +template +inline void AGEMOEA::SurvivalScoreAssignment( + const std::vector& front, + const arma::Col& idealPoint, + std::vector>& calculatedObjectives, + std::vector& survivalScore, + arma::Col& normalize, + double& dimension, + size_t fNum) +{ + typedef typename MatType::elem_type ElemType; + + // Calculations for the first front. + if (fNum == 0) + { + if (front.size() < numObjectives) + { + dimension = 1; + arma::Row extreme(numObjectives, arma::fill::zeros); + NormalizeFront(calculatedObjectives, normalize, front, extreme); + return; + } + + for (size_t index = 1; index < front.size(); index++) + { + calculatedObjectives[front[index]] = calculatedObjectives[front[index]] + - idealPoint; + } + + arma::Row extreme(numObjectives, arma::fill::zeros); + FindExtremePoints(extreme, calculatedObjectives, front); + NormalizeFront(calculatedObjectives, normalize, front, extreme); + + for (size_t index = 0; index < front.size(); index++) + { + calculatedObjectives[front[index]] = calculatedObjectives[front[index]] + / normalize; + } + + dimension = GetGeometry(calculatedObjectives, front, + extreme); + + std::set selected; + std::set remaining; + + // Create the selected and remaining sets. + for (size_t index: extreme) + { + selected.insert(index); + survivalScore[front[index]] = arma::datum::inf; + } + for (size_t i = 0; i < front.size(); i++) + { + if (selected.count(i) == 0) + { + remaining.insert(i); + } + } + + arma::Mat pairwise(front.size(), front.size(), arma::fill::zeros); + PairwiseDistance(pairwise,calculatedObjectives,front,dimension); + arma::Row proximity(front.size(), + arma::fill::zeros); + arma::Row diversity(front.size(), + arma::fill::zeros); + arma::Row value(front.size(), + arma::fill::zeros); + + // Calculate the diversity and proximity score. + for (size_t i = 0; i < front.size(); i++) + { + proximity[i] = std::pow(arma::accu(arma::pow( + arma::abs(calculatedObjectives[front[i]]), dimension)), 1.0 / dimension); + } + + while (remaining.size() > 0) + { + std::set::iterator it; + for (it = remaining.begin(); it != remaining.end(); it++) + { + diversity[*it] = DiversityScore(selected, pairwise, *it); + value[*it] = diversity[*it] / proximity[*it]; + } + size_t index = arma::index_max(value); + survivalScore[front[index]] = value[index]; + value[index] = - 1; + selected.insert(index); + remaining.erase(index); + } + } + + // Calculations for the other fronts. + else + { + for (size_t i = 0; i < front.size(); i++) + { + calculatedObjectives[front[i]] = calculatedObjectives[front[i]] / normalize; + survivalScore[front[i]] = std::pow(arma::accu(arma::pow(arma::abs( + calculatedObjectives[front[i]] - idealPoint), dimension)), + 1.0 / dimension); + } + + } +} + +//! Comparator for survival score based sorting. +template +inline bool AGEMOEA::SurvivalScoreOperator( + size_t idxP, + size_t idxQ, + const std::vector& ranks, + const std::vector& survivalScore) +{ + if (ranks[idxP] < ranks[idxQ]) + return true; + else if (ranks[idxP] == ranks[idxQ] && survivalScore[idxP] > survivalScore[idxQ]) + return true; + + return false; +} + +} // namespace ens + +#endif diff --git a/include/ensmallen_bits/utility/indicators/igd.hpp b/include/ensmallen_bits/utility/indicators/igd.hpp new file mode 100644 index 000000000..3adf1d772 --- /dev/null +++ b/include/ensmallen_bits/utility/indicators/igd.hpp @@ -0,0 +1,95 @@ +/** + * @file igd.hpp + * @author Satyam Shukla + * + * Inverse Generational Distance Plus (IGD) indicator. + * + * 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_INDICATORS_IGD_HPP +#define ENSMALLEN_INDICATORS_IGD_HPP + +namespace ens { + +/** + * The inverted generational distance( IGD) is a metric for assessing the quality + * of approximations to the Pareto front obtained by multi-objective optimization + * algorithms.The IGD indicator returns the average distance from each point in + * the reference front to the nearest point to it's solution. + * + * \f[ d(z,a) = \sqrt{\sum_{i = 1}^{n}(a_i - z_i)^2 \ } \ + * \f] + * + * For more information see: + * + * @code + * @inproceedings{coello2004study, + * title={A study of the parallelization of a coevolutionary multi-objective evolutionary algorithm}, + * author={Coello Coello, Carlos A and Reyes Sierra, Margarita}, + * booktitle={MICAI 2004: Advances in Artificial Intelligence: Third Mexican International Conference on Artificial Intelligence, Mexico City, Mexico, April 26-30, 2004. Proceedings 3}, + * pages={688--697}, + * year={2004}, + * organization={Springer} + * } + * @endcode + */ + class IGD + { + public: + /** + * Default constructor does nothing, but is required to satisfy the Indicator + * policy. + */ + IGD() { } + + /** + * Find the IGD value of the front with respect to the given reference + * front. + * + * @tparam CubeType The cube data type of front. + * @param front The given approximation front. + * @param referenceFront The given reference front. + * @param p The power constant in the distance formula. + * @return The IGD value of the front. + */ + template + static typename CubeType::elem_type Evaluate(const CubeType& front, + const CubeType& referenceFront, + double p) + { + // Convenience typedefs. + typedef typename CubeType::elem_type ElemType; + ElemType igd = 0; + for (size_t i = 0; i < referenceFront.n_slices; i++) + { + ElemType min = std::numeric_limits::max(); + for (size_t j = 0; j < front.n_slices; j++) + { + ElemType dist = 0; + for (size_t k = 0; k < front.slice(j).n_rows; k++) + { + ElemType z = referenceFront(k, 0, i); + ElemType a = front(k, 0, j); + // Assuming minimization of all objectives. + //! IGD does not clip negative differences to 0 + dist += std::pow(a - z, 2); + } + dist = std::sqrt(dist); + if (dist < min) + min = dist; + } + igd += std::pow(min,p); + } + igd /= referenceFront.n_slices; + igd = std::pow(igd, 1.0 / p); + return igd; + } + }; + +} // namespace ens + +#endif \ No newline at end of file diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index df5938844..515247ee0 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -24,12 +24,14 @@ set(ENSMALLEN_TESTS_SOURCES gradient_descent_test.cpp grid_search_test.cpp iqn_test.cpp + indicators_test.cpp katyusha_test.cpp lbfgs_test.cpp line_search_test.cpp lookahead_test.cpp lrsdp_test.cpp moead_test.cpp + agemoea_test.cpp momentum_sgd_test.cpp nesterov_momentum_sgd_test.cpp nsga2_test.cpp diff --git a/tests/agemoea_test.cpp b/tests/agemoea_test.cpp new file mode 100644 index 000000000..406975acb --- /dev/null +++ b/tests/agemoea_test.cpp @@ -0,0 +1,335 @@ +/** + * @file agemoea_test.cpp + * @author Satyam Shukla + * + * 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. + */ + +#include +#include "catch.hpp" +#include "test_function_tools.hpp" + +using namespace ens; +using namespace ens::test; +using namespace std; + +/** + * Checks if low <= value <= high. Used by MOEADFonsecaFlemingTest. + * + * @param value The value being checked. + * @param low The lower bound. + * @param high The upper bound. + * @param roundoff To round off precision. + * @tparam The type of elements in the population set. + * @return true if value lies in the range [low, high]. + * @return false if value does not lie in the range [low, high]. + */ +template +bool IsInBounds(const ElemType& value, + const ElemType& low, + const ElemType& high, + const ElemType& roundoff) +{ + return !(value < (low - roundoff)) && !((high + roundoff) < value); +} + +/** + * Optimize for the Schaffer N.1 function using AGE-MOEA optimizer. + * Tests for data of type double. + */ +TEST_CASE("AGEMOEASchafferN1DoubleTest", "[AGEMOEATest]") +{ + SchafferFunctionN1 SCH; + const double lowerBound = -1000; + const double upperBound = 1000; + const double expectedLowerBound = 0.0; + const double expectedUpperBound = 2.0; + + AGEMOEA opt(20, 500, 0.6, 20, 1e-6, 20, lowerBound, upperBound); + + typedef decltype(SCH.objectiveA) ObjectiveTypeA; + typedef decltype(SCH.objectiveB) ObjectiveTypeB; + + // We allow a few trials in case of poor convergence. + bool success = false; + for (size_t trial = 0; trial < 3; ++trial) + { + arma::mat coords = SCH.GetInitialPoint(); + std::tuple objectives = SCH.GetObjectives(); + + opt.Optimize(objectives, coords); + arma::cube paretoSet = opt.ParetoSet(); + + bool allInRange = true; + + for (size_t solutionIdx = 0; solutionIdx < paretoSet.n_slices; ++solutionIdx) + { + double val = arma::as_scalar(paretoSet.slice(solutionIdx)); + if (!IsInBounds(val, expectedLowerBound, expectedUpperBound, 0.1)) + { + allInRange = false; + break; + } + } + + if (allInRange) + { + success = true; + break; + } + } + + REQUIRE(success == true); +} + +/** + * Optimize for the Schaffer N.1 function using AGE-MOEA optimizer. + * Tests for data of type double. + */ +TEST_CASE("AGEMOEASchafferN1TestVectorDoubleBounds", "[AGEMOEATest]") +{ + // This test can be a little flaky, so we try it a few times. + SchafferFunctionN1 SCH; + const arma::vec lowerBound = {-1000}; + const arma::vec upperBound = {1000}; + const double expectedLowerBound = 0.0; + const double expectedUpperBound = 2.0; + + AGEMOEA opt(20, 500, 0.6, 20, 1e-6, 20, lowerBound, upperBound); + + typedef decltype(SCH.objectiveA) ObjectiveTypeA; + typedef decltype(SCH.objectiveB) ObjectiveTypeB; + + bool success = false; + for (size_t trial = 0; trial < 3; ++trial) + { + arma::mat coords = SCH.GetInitialPoint(); + std::tuple objectives = SCH.GetObjectives(); + + opt.Optimize(objectives, coords); + arma::cube paretoSet = opt.ParetoSet(); + + bool allInRange = true; + + for (size_t solutionIdx = 0; solutionIdx < paretoSet.n_slices; ++solutionIdx) + { + double val = arma::as_scalar(paretoSet.slice(solutionIdx)); + if (!IsInBounds(val, expectedLowerBound, expectedUpperBound, 0.1)) + { + allInRange = false; + break; + } + } + + if (allInRange) + { + success = true; + break; + } + } + + REQUIRE(success == true); +} + +/** + * Optimize for the Fonseca Fleming function using AGE-MOEA optimizer. + * Tests for data of type double. + */ +TEST_CASE("AGEMOEAFonsecaFlemingDoubleTest", "[AGEMOEATest]") +{ + FonsecaFlemingFunction FON; + const double lowerBound = -4; + const double upperBound = 4; + const double expectedLowerBound = -1.0 / sqrt(3); + const double expectedUpperBound = 1.0 / sqrt(3); + + AGEMOEA opt(20, 500, 0.6, 20, 1e-6, 20, lowerBound, upperBound); + typedef decltype(FON.objectiveA) ObjectiveTypeA; + typedef decltype(FON.objectiveB) ObjectiveTypeB; + + bool success = false; + for (size_t trial = 0; trial < 6; ++trial) + { + arma::mat coords = FON.GetInitialPoint(); + std::tuple objectives = FON.GetObjectives(); + + opt.Optimize(objectives, coords); + arma::cube paretoSet = opt.ParetoSet(); + + bool allInRange = true; + + for (size_t solutionIdx = 0; solutionIdx < paretoSet.n_slices; ++solutionIdx) + { + const arma::mat solution = paretoSet.slice(solutionIdx); + double valX = arma::as_scalar(solution(0)); + double valY = arma::as_scalar(solution(1)); + double valZ = arma::as_scalar(solution(2)); + + if (!IsInBounds(valX, expectedLowerBound, expectedUpperBound, 0.3) || + !IsInBounds(valY, expectedLowerBound, expectedUpperBound, 0.3) || + !IsInBounds(valZ, expectedLowerBound, expectedUpperBound, 0.3)) + { + allInRange = false; + break; + } + } + if(allInRange == true) + { + success = true; + break; + } + } + + REQUIRE(success == true); +} + +/** + * Optimize for the Fonseca Fleming function using AGE-MOEA optimizer. + * Tests for data of type float. + */ +TEST_CASE("AGEMOEAFonsecaFlemingTestVectorFloatBounds", "[AGEMOEATest]") +{ + FonsecaFlemingFunction FON; + const arma::vec lowerBound = {-4, -4, -4}; + const arma::vec upperBound = {4, 4, 4}; + const float expectedLowerBound = -1.0 / sqrt(3); + const float expectedUpperBound = 1.0 / sqrt(3); + + AGEMOEA opt(20, 300, 0.6, 20, 1e-6, 20, lowerBound, upperBound); + typedef decltype(FON.objectiveA) ObjectiveTypeA; + typedef decltype(FON.objectiveB) ObjectiveTypeB; + + bool success = false; + for (size_t trial = 0; trial < 6; ++trial) + { + arma::fmat coords = FON.GetInitialPoint(); + std::tuple objectives = FON.GetObjectives(); + + opt.Optimize(objectives, coords); + arma::fcube paretoSet = arma::conv_to::from(opt.ParetoSet()); + + bool allInRange = true; + for (size_t solutionIdx = 0; solutionIdx < paretoSet.n_slices; ++solutionIdx) + { + const arma::fmat solution = paretoSet.slice(solutionIdx); + float valX = arma::as_scalar(solution(0)); + float valY = arma::as_scalar(solution(1)); + float valZ = arma::as_scalar(solution(2)); + + if (!IsInBounds(valX, expectedLowerBound, expectedUpperBound, 0.3) || + !IsInBounds(valY, expectedLowerBound, expectedUpperBound, 0.3) || + !IsInBounds(valZ, expectedLowerBound, expectedUpperBound, 0.3)) + { + allInRange = false; + } + } + if (allInRange) + { + success = true; + break; + } + } + REQUIRE(success == true); +} + +TEST_CASE("AGEMOEAZDTONETest", "[AGEMOEATest]") +{ + //! Parameters taken from original ZDT Paper. + ZDT1<> ZDT_ONE(100); + const double lowerBound = 0; + const double upperBound = 1; + + AGEMOEA opt(20, 300, 0.6, 20, 1e-6, 20, lowerBound, upperBound); + + typedef decltype(ZDT_ONE.objectiveF1) ObjectiveTypeA; + typedef decltype(ZDT_ONE.objectiveF2) ObjectiveTypeB; + + const size_t trials = 8; + for (size_t trial = 0; trial < trials; ++trial) + { + arma::mat coords = ZDT_ONE.GetInitialPoint(); + std::tuple objectives = + ZDT_ONE.GetObjectives(); + + opt.Optimize(objectives, coords); + + //! Refer the ZDT_ONE implementation for g objective implementation. + //! The optimal g value is taken from the docs of ZDT_ONE. + size_t numVariables = coords.size(); + double sum = arma::accu(coords(arma::span(1, numVariables - 1), 0)); + const double g = 1.0 + 9.0 * sum / (static_cast(numVariables - 1)); + if (trial < trials - 1 && g != Approx(1.0).margin(0.99)) + continue; + + REQUIRE(g == Approx(1.0).margin(0.99)); + break; + } +} + +/** + * Check if the final population lies in the optimal region in variable space. + * + * @param paretoSet The final population in variable space. + */ +bool AVariableBoundsCheck(const arma::cube& paretoSet) +{ + const arma::mat regions{ + {0.0, 0.182228780, 0.4093136748, + 0.6183967944, 0.8233317983}, + {0.0830015349, 0.2577623634, 0.4538821041, + 0.6525117038, 0.8518328654} + }; + double notInBounds = 0; + for (size_t pointIdx = 0; pointIdx < paretoSet.n_slices; ++pointIdx) + { + const arma::mat& point = paretoSet.slice(pointIdx); + const double firstVariable = point(0, 0); + + const bool notInRegion0 = !IsInBounds(firstVariable, regions(0, 0), regions(1, 0), 3e-2); + const bool notInRegion1 = !IsInBounds(firstVariable, regions(0, 1), regions(1, 1), 3e-2); + const bool notInRegion2 = !IsInBounds(firstVariable, regions(0, 2), regions(1, 2), 3e-2); + const bool notInRegion3 = !IsInBounds(firstVariable, regions(0, 3), regions(1, 3), 3e-2); + const bool notInRegion4 = !IsInBounds(firstVariable, regions(0, 4), regions(1, 4), 3e-2); + + if (notInRegion0 && notInRegion1 && notInRegion2 && notInRegion3 && notInRegion4) + { + notInBounds++; + } + } + + notInBounds = notInBounds / paretoSet.n_slices; + return notInBounds < 0.80; +} + +/** + * Test AGEMOEA against the third problem of ZDT Test Suite. ZDT-3 is a 30 + * variable-2 objective problem with disconnected Pareto Fronts. + */ +TEST_CASE("AGEMOEADIRICHLETZDT3Test", "[AGEMOEADTest]") +{ + //! Parameters taken from original ZDT Paper. + ZDT3<> ZDT_THREE(300); + const double lowerBound = 0; + const double upperBound = 1; + + AGEMOEA opt(50, 500, 0.8, 20, 1e-6, 20, lowerBound, upperBound); + + typedef decltype(ZDT_THREE.objectiveF1) ObjectiveTypeA; + typedef decltype(ZDT_THREE.objectiveF2) ObjectiveTypeB; + bool success = true; + for (size_t tries = 0; tries < 2; tries++) + { + arma::mat coords = ZDT_THREE.GetInitialPoint(); + std::tuple objectives = ZDT_THREE.GetObjectives(); + + opt.Optimize(objectives, coords); + + const arma::cube& finalPopulation = opt.ParetoSet(); + success = AVariableBoundsCheck(finalPopulation); + if (success){ break;} + } + REQUIRE(success); +} \ No newline at end of file diff --git a/tests/indicators_test.cpp b/tests/indicators_test.cpp index 3f641fcbe..478a1ea2c 100644 --- a/tests/indicators_test.cpp +++ b/tests/indicators_test.cpp @@ -2,7 +2,7 @@ * @file indicators_test.cpp * @author Nanubala Gnana Sai * - * Test file for all the indicators: Epsilon, IGD+. + * Test file for all the indicators: Epsilon, IGD+, IGD. * * 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 @@ -17,7 +17,7 @@ using namespace ens; using namespace ens::test; /** - * Calculates the Epsilon metric for the pair of fronts. + * Calculates the Epsilon performance indicator for the pair of fronts. * Tests for data of type double. * The reference numerical results have been taken from hand calculated values. * Refer the IPynb notebook in https://github.com/mlpack/ensmallen/pull/285 @@ -37,7 +37,7 @@ TEST_CASE("EpsilonDoubleTest", "[IndicatorsTest]") } /** - * Calculates the Epsilon metric for the pair of fronts. + * Calculates the Epsilon performance indicator for the pair of fronts. * Tests for data of type float. * The reference numerical results have been taken from hand calculated values. * Refer the IPynb notebook in https://github.com/mlpack/ensmallen/pull/285 @@ -57,10 +57,49 @@ TEST_CASE("EpsilonFloatTest", "[IndicatorsTest]") } /** - * Calculates the IGD+ metric for the pair of fronts. + * Calculates the IGD performance indicator for the pair of fronts. * Tests for data of type double. * The reference numerical results have been taken from hand calculated values. - * Refer the IPynb notebook in https://github.com/mlpack/ensmallen/pull/285 + * Refer the IPynb notebook in https://github.com/mlpack/ensmallen/pull/ + * for more. + */ +TEST_CASE("IGDDoubleTest", "[IndicatorsTest]") +{ + arma::cube referenceFront(2, 1, 3); + double tol = 1e-10; + referenceFront.slice(0) = arma::vec{0.11111111, 0.75039705}; + referenceFront.slice(1) = arma::vec{0.22222222, 0.60558677}; + referenceFront.slice(2) = arma::vec{0.33333333, 0.49446993}; + arma::cube front = referenceFront * 1.1; + double igd = IGD::Evaluate(front, referenceFront, 1); + REQUIRE(igd == Approx(0.06666608265617857).margin(tol)); +} + +/** + * Calculates the IGD performance indicator for the pair of fronts. + * Tests for data of type float. + * The reference numerical results have been taken from hand calculated values. + * Refer the IPynb notebook in https://github.com/mlpack/ensmallen/pull/399 + * for more. + */ +TEST_CASE("IGDFloatTest", "[IndicatorsTest]") +{ + arma::fcube referenceFront(2, 1, 3); + float tol = 1e-10; + referenceFront.slice(0) = arma::fvec{0.11111111f, 0.75039705f}; + referenceFront.slice(1) = arma::fvec{0.22222222f, 0.60558677f}; + referenceFront.slice(2) = arma::fvec{0.33333333f, 0.49446993f}; + arma::fcube front = referenceFront * 1.1; + float igd = IGD::Evaluate(front, referenceFront, 1); + + REQUIRE(igd == Approx(0.06666608265617857).margin(tol)); +} + +/** + * Calculates the IGD+ performance indicator for the pair of fronts. + * Tests for data of type double. + * The reference numerical results have been taken from hand calculated values. + * Refer the IPynb notebook in https://github.com/mlpack/ensmallen/pull/399 * for more. */ TEST_CASE("IGDPlusDoubleTest", "[IndicatorsTest]") @@ -77,7 +116,7 @@ TEST_CASE("IGDPlusDoubleTest", "[IndicatorsTest]") } /** - * Calculates the IGD+ metric for the pair of fronts. + * Calculates the IGD+ performance indicator for the pair of fronts. * Tests for data of type float. * The reference numerical results have been taken from hand calculated values. * Refer the IPynb notebook in https://github.com/mlpack/ensmallen/pull/285