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

using System;
using System.Collections.Generic;
using System.Linq;
using HeuristicLab.Core;
using HeuristicLab.Optimization;

namespace HeuristicLab.Algorithms.NSGA3
{
    public static class NSGA3Selection
    {
        private const double EPSILON = 10e-6; // a tiny number that is greater than 0

        /// <summary>
        /// TODO: Description If the fitnesses of all solutions in <paramref name="rt" /> do not
        /// have the same number of objectives, undefined behavior happens.
        /// </summary>
        /// <param name="rt"></param>
        /// <param name="referencePoints"></param>
        /// <param name="maximization"></param>
        /// <returns></returns>
        public static List<Solution> SelectSolutionsForNextGeneration(List<Solution> rt, List<ReferencePoint> referencePoints, bool[] maximization, IRandom random)
        {
            int N = rt.Count / 2; // number of solutions in a generation
            int numberOfObjectives = rt.First().Fitness.Length;
            List<Solution> st = new List<Solution>();

            // Do non-dominated sort
            var qualities = Utility.ToFitnessMatrix(rt);
            // compute the pareto fronts using the DominationCalculator and discard the qualities
            // part in the inner tuples
            var fronts = DominationCalculator<Solution>.CalculateAllParetoFronts(rt.ToArray(), qualities, maximization, out int[] rank, true)
                .Select(list => new List<Solution>(list.Select(pair => pair.Item1))).ToList();

            int lf = -1; // last front to be included
            for (int i = 0; i < fronts.Count && st.Count < N; i++)
            {
                st = Utility.Concat(st, fronts[i]);
                lf = i;
            }
            // remove useless fronts
            fronts.RemoveRange(lf + 1, fronts.Count - lf - 1);

            List<Solution> nextGeneration = new List<Solution>();
            if (st.Count == N) // 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 < lf; f++)
                    nextGeneration = Utility.Concat(nextGeneration, fronts[f]);
                int k = N - nextGeneration.Count;
                Normalize(st, numberOfObjectives);
                Associate(fronts, referencePoints);
                List<Solution> solutionsToAdd = Niching(k, referencePoints, random);
                nextGeneration = Utility.Concat(nextGeneration, solutionsToAdd);
            }

            return nextGeneration;
        }

        private static void Normalize(List<Solution> st, int numberOfObjectives)
        {
            // Find the ideal point
            double[] idealPoint = new double[numberOfObjectives];
            for (int j = 0; j < numberOfObjectives; j++)
            {
                // Compute ideal point
                idealPoint[j] = Utility.Min(s => s.Fitness[j], st);

                // Translate objectives
                foreach (var solution in st)
                    solution.Fitness[j] -= idealPoint[j];
            }

            // Find the extreme points
            Solution[] extremePoints = new 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;

                extremePoints[j] = Utility.ArgMin(s => ASF(s.Fitness, weights), st);
            }

            // Compute intercepts
            List<double> intercepts = GetIntercepts(extremePoints.ToList());

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

        private static void NormalizeObjectives(List<Solution> st, List<double> intercepts, double[] idealPoint)
        {
            int numberOfObjectives = st.First().Fitness.Length;

            foreach (var solution in st)
            {
                for (int i = 0; i < numberOfObjectives; i++)
                {
                    if (Math.Abs(intercepts[i] - idealPoint[i]) > EPSILON)
                    {
                        solution.Fitness[i] = solution.Fitness[i] / (intercepts[i] - idealPoint[i]);
                    }
                    else
                    {
                        solution.Fitness[i] = solution.Fitness[i] / EPSILON;
                    }
                }
            }
        }

        private static 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.Fitness), referencePoints);
                    // associated reference point
                    var arp = rpAndDist.Item1;
                    // distance to that reference point
                    var dist = rpAndDist.Item2;

                    if (f + 1 != fronts.Count)
                    {
                        // Todo: Add member for reference point on index min_rp
                        arp.NumberOfAssociatedSolutions++;
                    }
                    else
                    {
                        // Todo: Add potential member for reference point on index min_rp
                        arp.AddPotentialAssociatedSolution(solution, dist);
                    }
                }
            }
        }

        private static 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 static 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);

            // the reference points that share the number of associated solutions where that number
            // is equal to minNumber
            List<ReferencePoint> minAssociatedReferencePoints = new List<ReferencePoint>();
            foreach (var referencePoint in referencePoints)
                if (referencePoint.NumberOfAssociatedSolutions == minNumber)
                    minAssociatedReferencePoints.Add(referencePoint);

            if (minAssociatedReferencePoints.Count > 1)
                return minAssociatedReferencePoints[random.Next(minAssociatedReferencePoints.Count)];
            else
                return minAssociatedReferencePoints.Single();
        }

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

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

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

        private static double ASF(double[] x, double[] weight)
        {
            int numberOfObjectives = x.Length;
            List<int> dimensions = new List<int>();
            for (int i = 0; i < numberOfObjectives; i++)
                dimensions.Add(i);

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

        private static List<double> GetIntercepts(List<Solution> extremePoints)
        {
            int numberOfObjectives = extremePoints.First().Fitness.Length;

            // 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++)
                {
                    // maybe todo: override Equals method of solution?
                    duplicate = extremePoints[i].Equals(extremePoints[j]);
                }
            }

            List<double> intercepts = new List<double>();

            if (duplicate)
            { // cannot construct the unique hyperplane (this is a casual method to deal with the condition)
                for (int f = 0; f < numberOfObjectives; f++)
                {
                    // extreme_points[f] stands for the individual with the largest value of
                    // objective f
                    intercepts.Add(extremePoints[f].Fitness[f]);
                }
            }
            else
            {
                // Find the equation of the hyperplane
                List<double> b = new List<double>(); //(pop[0].objs().size(), 1.0);
                for (int i = 0; i < numberOfObjectives; i++)
                {
                    b.Add(1.0);
                }

                List<List<double>> a = new List<List<double>>();
                foreach (Solution s in extremePoints)
                {
                    List<double> aux = new List<double>();
                    for (int i = 0; i < numberOfObjectives; i++)
                        aux.Add(s.Fitness[i]);
                    a.Add(aux);
                }
                List<double> x = GaussianElimination(a, b);

                // Find intercepts
                for (int f = 0; f < numberOfObjectives; f++)
                {
                    intercepts.Add(1.0 / x[f]);
                }
            }

            return intercepts;
        }

        private static 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;
        }
    }
}