source: branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/NSGA3Selection.cs

using System;
using System.Collections.Generic;
using System.Linq;
using HEAL.Attic;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Optimization;

namespace HeuristicLab.Algorithms.NSGA3
{
    [StorableType("53158EA5-4D1C-42DB-90F5-EC555EC0A450")]
    public class NSGA3Selection : IDeepCloneable
    {
        private const double EPSILON = 10e-6; // a tiny number that is greater than 0

        [Storable]
        private int numberOfObjectives;

        [Storable]
        private bool[] maximization;

        [Storable]
        private int populationSize;

        [Storable]
        private IRandom random;

        [Storable]
        private bool dominateOnEqualQualities;

        [StorableConstructor]
        public NSGA3Selection(StorableConstructorFlag _) { }

        public NSGA3Selection(int numberOfObjectives, bool[] maximization, int populationSize, IRandom random, bool dominateOnEqualQualities)
        {
            if (maximization == null) throw new ArgumentNullException(nameof(maximization));
            if (random == null) throw new ArgumentNullException(nameof(random));
            if (numberOfObjectives != maximization.Length) throw new ArgumentException($"{nameof(numberOfObjectives)} and the length of {nameof(maximization)} do not match.");
            if (populationSize < 1) throw new ArgumentOutOfRangeException($"{nameof(populationSize)} has to be >= 1.");

            this.numberOfObjectives = numberOfObjectives;
            this.maximization = maximization;
            this.populationSize = populationSize;
            this.random = random;
            this.dominateOnEqualQualities = dominateOnEqualQualities;
        }

        public NSGA3Selection(NSGA3Selection other, Cloner cloner)
        {
            numberOfObjectives = other.numberOfObjectives;
            maximization = new bool[other.maximization.Length];
            Array.Copy(other.maximization, maximization, maximization.Length);
            populationSize = other.populationSize;
            random = cloner.Clone(other.random);
            dominateOnEqualQualities = other.dominateOnEqualQualities;
        }

        public IDeepCloneable Clone(Cloner cloner)
        {
            return new NSGA3Selection(this, cloner);
        }

        public object Clone()
        {
            return new Cloner().Clone(this);
        }

        /// <summary>
        /// Select the solutions for the next generation from <see cref="rt" /> according to the
        /// NSGA-III algorithm.
        /// </summary>
        /// <param name="rt">The current generation and its mutated children.</param>
        /// <param name="referencePoints"></param>
        /// <returns></returns>
        public List<Solution> SelectSolutionsForNextGeneration(List<Solution> rt, List<ReferencePoint> referencePoints)
        {
            List<Solution> st = new List<Solution>();

            // Do non-dominated sort
            var qualities = Utility.ToFitnessMatrix(rt);

            // compute the pareto fronts using the DominationCalculator
            List<List<Solution>> fronts = DominationCalculator<Solution>.CalculateAllParetoFronts(rt.ToArray(), qualities, maximization, out int[] rank, dominateOnEqualQualities)
                .Select(list => new List<Solution>(list.Select(pair => pair.Item1))).ToList();

            int last = 0; // last front to be included
            while (st.Count < populationSize)
            {
                st = Utility.Concat(st, fronts[last++]);
            }

            // remove useless fronts
            fronts.RemoveRange(last, fronts.Count - last);

            List<Solution> nextGeneration = new List<Solution>();
            if (st.Count == populationSize)
            { // no special selection needs to be done
                nextGeneration = st;
            }
            else
            {
                // Add every front to the next generation list except for the one that is added only partially
                nextGeneration = new List<Solution>();
                for (int f = 0; f < last - 1; f++)
                    nextGeneration = Utility.Concat(nextGeneration, fronts[f]);

                int k = populationSize - nextGeneration.Count;
                Normalize(rt, fronts);
                Associate(fronts, referencePoints);
                List<Solution> solutionsToAdd = Niching(k, referencePoints, random);
                nextGeneration = Utility.Concat(nextGeneration, solutionsToAdd);
            }

            return nextGeneration;
        }

        private void Normalize(List<Solution> solutions, List<List<Solution>> fronts)
        {
            // Find the ideal point
            double[] idealPoint = TranslateObjectives(solutions, fronts[0]);

            // Find the extreme points
            List<Solution> extremePoints = GetExtremePoints(fronts[0]);

            // Compute intercepts
            double[] intercepts = GetIntercepts(extremePoints, fronts[0]);

            // Normalize objectives
            NormalizeObjectives(solutions, intercepts, idealPoint);
        }

        private double[] TranslateObjectives(List<Solution> solutions, List<Solution> firstFront)
        {
            double[] idealPoint = new double[numberOfObjectives];
            foreach (var solution in solutions)
            {
                if (solution.NormalizedObjectives == null) solution.NormalizedObjectives = new double[numberOfObjectives];
            }

            for (int j = 0; j < numberOfObjectives; j++)
            {
                // Compute ideal point
                idealPoint[j] = Utility.Min(s => s.Objectives[j], firstFront);

                // Translate objectives and save the modified fitness values in NormalizedObjectives
                foreach (var solution in solutions)
                    solution.NormalizedObjectives[j] = solution.Objectives[j] - idealPoint[j];
            }

            return idealPoint;
        }

        private List<Solution> GetExtremePoints(List<Solution> firstFront)
        {
            List<Solution> extremePoints = new List<Solution>(numberOfObjectives);
            for (int j = 0; j < numberOfObjectives; j++)
            {
                // Compute extreme points
                double[] weights = new double[numberOfObjectives];
                for (int i = 0; i < numberOfObjectives; i++) weights[i] = EPSILON;
                weights[j] = 1;

                var extremePoint = Utility.ArgMin(s => ASF(s.NormalizedObjectives, weights), firstFront);
                extremePoints.Add(extremePoint);
            }

            return extremePoints;
        }

        private double ASF(double[] convertedFitness, double[] weight)
        {
            var dims = Enumerable.Range(0, numberOfObjectives);

            return Utility.Max((dim) => convertedFitness[dim] / weight[dim], dims);
        }

        private double[] GetIntercepts(List<Solution> extremePoints, List<Solution> solutions)
        {
            // Check whether there are duplicate extreme points. This might happen but the original
            // paper does not mention how to deal with it.
            bool duplicate = false;
            for (int i = 0; !duplicate && i < extremePoints.Count; i++)
            {
                for (int j = i + 1; !duplicate && j < extremePoints.Count; j++)
                {
                    duplicate = extremePoints[i] == extremePoints[j];
                }
            }

            double[] intercepts = new double[numberOfObjectives];

            bool negative_intercept = false;
            if (!duplicate)
            {
                // Find the equation of the hyperplane
                List<double> b = new List<double>(Enumerable.Repeat(1.0, numberOfObjectives));

                List<List<double>> a = new List<List<double>>(extremePoints.Count);

                for (int i = 0; i < extremePoints.Count; i++)
                {
                    List<double> convertedFitness = new List<double>(extremePoints[i].NormalizedObjectives);
                    a.Add(convertedFitness);
                }
                List<double> x = GaussianElimination(a, b);

                // Find intercepts
                for (int i = 0; i < intercepts.Length; i++)
                {
                    intercepts[i] = 1.0 / x[i];

                    if (x[i] < 0)
                    {
                        negative_intercept = true;
                        break;
                    }
                }
            }

            // follow the method in Yuan et al. (GECCO 2015)
            if (duplicate || negative_intercept)
            {
                double[] maxObjectives = FindMaxObjectives(solutions);
                for (int i = 0; i < intercepts.Length; i++)
                    intercepts[i] = maxObjectives[i];
            }

            return intercepts;
        }

        private List<double> GaussianElimination(List<List<double>> a, List<double> b)
        {
            List<double> x = new List<double>();

            int n = a.Count;
            for (int i = 0; i < n; i++)
                a[i].Add(b[i]);

            for (int @base = 0; @base < n - 1; @base++)
                for (int target = @base + 1; target < n; target++)
                {
                    double ratio = a[target][@base] / a[@base][@base];
                    for (int term = 0; term < a[@base].Count; term++)
                        a[target][term] = a[target][term] - a[@base][term] * ratio;
                }

            for (int i = 0; i < n; i++)
                x.Add(0.0);

            for (int i = n - 1; i >= 0; i--)
            {
                for (int known = i + 1; known < n; known++)
                    a[i][n] = a[i][n] - a[i][known] * x[known];
                x[i] = a[i][n] / a[i][i];
            }

            return x;
        }

        private double[] FindMaxObjectives(List<Solution> solutions)
        {
            double[] maxPoint = new double[numberOfObjectives];
            for (int i = 0; i < maxPoint.Length; i++) maxPoint[i] = double.MinValue;

            for (int f = 0; f < numberOfObjectives; f++)
            {
                maxPoint[f] = Utility.Max(s => s.Objectives[f], solutions);
            }

            return maxPoint;
        }

        private void NormalizeObjectives(List<Solution> solutions, double[] intercepts, double[] idealPoint)
        {
            foreach (var solution in solutions)
            {
                for (int i = 0; i < numberOfObjectives; i++)
                {
                    if (Math.Abs(intercepts[i] - idealPoint[i]) > 10e-10)
                    {
                        solution.NormalizedObjectives[i] = solution.NormalizedObjectives[i] / (intercepts[i] - idealPoint[i]);
                    }
                    else
                    {
                        solution.NormalizedObjectives[i] = solution.NormalizedObjectives[i] / 10e-10;
                    }
                }
            }
        }

        private void Associate(List<List<Solution>> fronts, List<ReferencePoint> referencePoints)
        {
            for (int f = 0; f < fronts.Count; f++)
            {
                foreach (var solution in fronts[f])
                {
                    // find reference point for which the perpendicular distance to the current
                    // solution is the lowest
                    var rpAndDist = Utility.MinArgMin(rp => GetPerpendicularDistance(rp.Values, solution.NormalizedObjectives), referencePoints);
                    // associated reference point
                    var arp = rpAndDist.Item1;
                    // distance to that reference point
                    var dist = rpAndDist.Item2;

                    if (f + 1 != fronts.Count)
                    {
                        arp.NumberOfAssociatedSolutions++;
                    }
                    else
                    {
                        arp.AddPotentialAssociatedSolution(solution, dist);
                    }
                }
            }
        }

        private double GetPerpendicularDistance(double[] direction, double[] point)
        {
            double numerator = 0;
            double denominator = 0;
            for (int i = 0; i < direction.Length; i++)
            {
                numerator += direction[i] * point[i];
                denominator += Math.Pow(direction[i], 2);
            }
            double k = numerator / denominator;

            double d = 0;
            for (int i = 0; i < direction.Length; i++)
            {
                d += Math.Pow(k * direction[i] - point[i], 2);
            }
            return Math.Sqrt(d);
        }

        private List<Solution> Niching(int k, List<ReferencePoint> referencePoints, IRandom random)
        {
            List<Solution> solutions = new List<Solution>();
            while (solutions.Count != k)
            {
                ReferencePoint min_rp = FindNicheReferencePoint(referencePoints, random);

                Solution chosen = SelectClusterMember(min_rp);
                if (chosen == null)
                {
                    referencePoints.Remove(min_rp);
                }
                else
                {
                    min_rp.NumberOfAssociatedSolutions++;
                    min_rp.RemovePotentialAssociatedSolution(chosen);
                    solutions.Add(chosen);
                }
            }

            return solutions;
        }

        private ReferencePoint FindNicheReferencePoint(List<ReferencePoint> referencePoints, IRandom random)
        {
            // the minimum number of associated solutions for a reference point over all reference points
            int minNumber = Utility.Min(rp => rp.NumberOfAssociatedSolutions, referencePoints);

            var minReferencePoints = new List<ReferencePoint>(referencePoints.Where(r => r.NumberOfAssociatedSolutions == minNumber));
            if (minReferencePoints.Count == 1) return minReferencePoints.Single();
            else return minReferencePoints[random.Next(minReferencePoints.Count)];
        }

        private Solution SelectClusterMember(ReferencePoint referencePoint)
        {
            Solution chosen = null;
            if (referencePoint.HasPotentialMember())
            {
                if (referencePoint.NumberOfAssociatedSolutions == 0)
                    chosen = referencePoint.FindClosestMember();
                else
                    chosen = referencePoint.RandomMember();
            }
            return chosen;
        }
    }
}