Changeset 17616

Timestamp:

06/21/20 16:34:08 (4 years ago)

Author:

dleko

Message:

#2825 Implement the first part of NSGA3.
Solutions are saved in a List.
Minor changes in Utility.cs.

Location:

branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3

Files:

• NSGA3.cs (modified) (6 diffs)
• Utility.cs (modified) (3 diffs)

branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/NSGA3.cs

@@ -25 +25 @@
     public class NSGA3 : BasicAlgorithm
     {
-        public override bool SupportsPause => false; // todo: make true
+        private const double EPSILON = 10e-6; // a tiny number that is greater than 0
+
+        public override bool SupportsPause => false;
 
         #region ProblemProperties
@@ -48 +50 @@
 
         [Storable]
-        private Solution[] solutions;
+        private List<Solution> solutions;
+
+        [Storable]
+        private List<List<Solution>> fronts;
 
         [Storable]
@@ -170 +175 @@
             // todo: don't forget to clone storable fields
             random = cloner.Clone(original.random);
-            solutions = original.solutions?.Select(cloner.Clone).ToArray();
+            solutions = new List<Solution>(original.solutions?.Select(cloner.Clone));
         }
 
@@ -192 +197 @@
         protected override void Run(CancellationToken cancellationToken)
         {
-            throw new NotImplementedException();
+            referencePoints = new List<ReferencePoint>(); // todo: use existing list of reference points
+            for (int iteration = 0; iteration < MaximumGenerations.Value; iteration++)
+                ToNextGeneration();
         }
 
@@ -207 +214 @@
         private void InitSolutions()
         {
-            int minBound = 0; // todo: find min inside Problem.Encoding
-            int maxBound = 1; // todo: find max inside Problem.Encoding
+            int minBound = 0;
+            int maxBound = 1;
 
             // Initialise solutions
-            solutions = new Solution[PopulationSize.Value];
+            solutions = new List<Solution>(PopulationSize.Value);
             for (int i = 0; i < PopulationSize.Value; i++)
             {
                 RealVector randomRealVector = new RealVector(Problem.Encoding.Length, random, minBound, maxBound);
 
-                solutions[i] = new Solution(StorableConstructorFlag.Default)
+                solutions.Add(new Solution(StorableConstructorFlag.Default)
                 {
                     Chromosome = randomRealVector
-                };
+                });
                 solutions[i].Fitness = Evaluate(solutions[i].Chromosome);
             }
@@ -260 +267 @@
         }
 
+        private List<Solution> ToNextGeneration()
+        {
+            List<Solution> st = new List<Solution>();
+            List<Solution> qt = Mutate(Recombine(solutions));
+            List<Solution> rt = Utility.Concat(solutions, qt);
+            List<Solution> nextGeneration;
+
+            // Do non-dominated sort
+            var qualities = Utility.ToFitnessMatrix(solutions);
+            // compute the pareto fronts using the DominationCalculator and discard the qualities
+            // part in the inner tuples
+            fronts = DominationCalculator<Solution>.CalculateAllParetoFronts(rt.ToArray(), qualities, Problem.Maximization, out int[] rank, true)
+                .Select(list => new List<Solution>(list.Select(pair => pair.Item1))).ToList();
+
+            int i = 0;
+            List<Solution> lf = null; // last front to be included
+            while (i < fronts.Count && st.Count < PopulationSize.Value)
+            {
+                lf = fronts[i];
+                st = Utility.Concat(st, lf);
+                i++;
+            }
+
+            if (st.Count == PopulationSize.Value) // no selection needs to be done
+                nextGeneration = st;
+            else
+            {
+                int l = i - 1;
+                nextGeneration = new List<Solution>();
+                for (int f = 0; f < l; f++)
+                    nextGeneration = Utility.Concat(nextGeneration, fronts[f]);
+                Normalize(st);
+                Associate();
+                throw new NotImplementedException();
+            }
+
+            throw new NotImplementedException();
+        }
+
+        private void Normalize(List<Solution> population)
+        {
+            // Find the ideal point
+            double[] idealPoint = new double[Problem.Encoding.Length];
+            for (int j = 0; j < Problem.Encoding.Length; j++)
+            {
+                // Compute ideal point
+                idealPoint[j] = Utility.Min(s => s.Fitness[j], solutions);
+
+                // Translate objectives
+                foreach (var solution in solutions)
+                    solution.Fitness[j] -= idealPoint[j];
+            }
+
+            // Find the extreme points
+            Solution[] extremePoints = new Solution[Problem.Encoding.Length];
+            for (int j = 0; j < Problem.Encoding.Length; j++)
+            {
+                // Compute extreme points
+                double[] weights = new double[Problem.Encoding.Length];
+                for (int i = 0; i < Problem.Encoding.Length; i++) weights[i] = EPSILON;
+                weights[j] = 1;
+                Func<Solution, double> func = s => ASF(s.Fitness, weights);
+                extremePoints[j] = Utility.ArgMin(func, solutions);
+            }
+
+            // Compute intercepts
+            List<double> intercepts = GetIntercepts(extremePoints.ToList());
+
+            // Normalize objectives
+            NormalizeObjectives(intercepts, idealPoint);
+
+            // Associate reference points to solutions
+            Associate();
+        }
+
+        private void NormalizeObjectives(List<double> intercepts, double[] idealPoint)
+        {
+            for (int f = 0; f < fronts.Count; f++)
+            {
+                foreach (var solution in fronts[f])
+                {
+                    for (int i = 0; i < Problem.Encoding.Length; 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 void Associate()
+        {
+            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);
+
+                    //// todo: use ArgMin here
+                    //int min_rp = -1;
+                    //double min_dist = double.MaxValue;
+                    //for (int r = 0; r < referencePoints.Count; r++)
+                    //{
+                    //    double d = GetPerpendicularDistance(referencePoints[r].Values, solution.Fitness);
+                    //    if (d < min_dist)
+                    //    {
+                    //        min_dist = d;
+                    //        min_rp = r;
+                    //    }
+                    //}
+                    if (f + 1 != fronts.Count)
+                    {
+                        // Todo: Add member for reference point on index min_rp
+                        throw new NotImplementedException();
+                    }
+                    else
+                    {
+                        // Todo: Add potential member for reference point on index min_rp
+                        throw new NotImplementedException();
+                    }
+                }
+            }
+        }
+
+        private 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 double ASF(double[] x, double[] weight)
+        {
+            List<int> dimensions = new List<int>();
+            for (int i = 0; i < Problem.Encoding.Length; i++) dimensions.Add(i);
+            Func<int, double> f = dim => x[dim] / weight[dim];
+            return Utility.Max(f, dimensions);
+        }
+
+        private List<double> GetIntercepts(List<Solution> extremePoints)
+        {
+            // 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 < Problem.Encoding.Length; 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 < Problem.Encoding.Length; 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 < Problem.Encoding.Length; i++)
+                        aux.Add(s.Fitness[i]);
+                    a.Add(aux);
+                }
+                List<double> x = GaussianElimination(a, b);
+
+                // Find intercepts
+                for (int f = 0; f < Problem.Encoding.Length; f++)
+                {
+                    intercepts.Add(1.0 / x[f]);
+                }
+            }
+
+            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 List<Solution> Recombine(List<Solution> solutions)
+        {
+            throw new NotImplementedException();
+        }
+
+        private List<Solution> Mutate(List<Solution> solutions)
+        {
+            throw new NotImplementedException();
+        }
+
         #endregion Private Methods
     }

branches/2825-NSGA3/HeuristicLab.Algorithms.NSGA3/3.3/Utility.cs

@@ -8 +8 @@
     internal static class Utility
     {
+        internal static List<T> Concat<T>(List<T> list1, List<T> list2)
+        {
+            List<T> resultList = new List<T>(list1.Count + list2.Count);
+
+            resultList.AddRange(list1);
+            resultList.AddRange(list2);
+
+            return resultList;
+        }
+
         /// <summary>
         /// Returns the number of possible combinations of size <paramref name="r" /> from a set of
@@ -42 +52 @@
         }
 
-        internal static double[,] ToArray(List<ReferencePoint> referencePoints)
+        internal static double[][] ToFitnessMatrix(this List<Solution> solutions)
+        {
+            double[][] data = new double[solutions.Count][];
+            for (int i = 0; i < solutions.Count; i++)
+            {
+                Solution solution = solutions[i];
+                data[i] = new double[solution.Fitness.Length];
+                for (int j = 0; j < solution.Fitness.Length; j++)
+                {
+                    data[i][j] = solution.Fitness[j];
+                }
+            }
+
+            return data;
+        }
+
+        internal static DoubleMatrix ToMatrix(this IEnumerable<IReadOnlyList<double>> data)
+        {
+            var d2 = data.ToArray();
+            var mat = new DoubleMatrix(d2.Length, d2[0].Count);
+            for (var i = 0; i < mat.Rows; i++)
+                for (var j = 0; j < mat.Columns; j++)
+                    mat[i, j] = d2[i][j];
+            return mat;
+        }
+
+        private static double[,] ToArray(List<ReferencePoint> referencePoints)
         {
             if (referencePoints == null || !referencePoints.Any()) throw new ArgumentException($"{nameof(referencePoints)} is null or empty");
@@ -56 +92 @@
         }
 
-        internal static DoubleMatrix ToMatrix(this IEnumerable<IReadOnlyList<double>> data)
+        internal static TOut Min<TIn, TOut>(Func<TIn, TOut> func, IEnumerable<TIn> inputs) where TOut : IComparable
         {
-            var d2 = data.ToArray();
-            var mat = new DoubleMatrix(d2.Length, d2[0].Count);
-            for (var i = 0; i < mat.Rows; i++)
-                for (var j = 0; j < mat.Columns; j++)
-                    mat[i, j] = d2[i][j];
-            return mat;
+            //var it = inputs.GetEnumerator();
+            //it.MoveNext();
+            //var minValue = func(it.Current);
+            //while (it.MoveNext())
+            //{
+            //    var currentValue = func(it.Current);
+            //    if (minValue.CompareTo(currentValue) > 0) minValue = currentValue;
+            //}
+
+            //return minValue;
+            return MinArgMin(func, inputs).Item2;
+        }
+
+        internal static TIn ArgMin<TIn, TOut>(Func<TIn, TOut> func, IEnumerable<TIn> inputs) where TOut : IComparable
+        {
+            return MinArgMin(func, inputs).Item1;
+        }
+
+        internal static Tuple<TIn, TOut> MinArgMin<TIn, TOut>(Func<TIn, TOut> func, IEnumerable<TIn> inputs) where TOut : IComparable
+        {
+            var it = inputs.GetEnumerator();
+            it.MoveNext();
+            var minArg = it.Current;
+            var minValue = func(it.Current);
+            while (it.MoveNext())
+            {
+                var currentArg = it.Current;
+                var currentValue = func(it.Current);
+                if (minValue.CompareTo(currentValue) > 0)
+                {
+                    minArg = currentArg;
+                    minValue = currentValue;
+                }
+            }
+
+            return Tuple.Create(minArg, minValue);
+        }
+
+        internal static TOut Max<TIn, TOut>(Func<TIn, TOut> func, IEnumerable<TIn> inputs) where TOut : IComparable
+        {
+            //var it = inputs.GetEnumerator();
+            //it.MoveNext();
+            //var maxValue = func(it.Current);
+            //while (it.MoveNext())
+            //{
+            //    var currentValue = func(it.Current);
+            //    if (maxValue.CompareTo(currentValue) < 0) maxValue = currentValue;
+            //}
+
+            //return maxValue;
+            return MaxArgMax(func, inputs).Item2;
+        }
+
+        internal static TIn ArgMax<TIn, TOut>(Func<TIn, TOut> func, IEnumerable<TIn> inputs) where TOut : IComparable
+        {
+            return MaxArgMax(func, inputs).Item1;
+        }
+
+        internal static Tuple<TIn, TOut> MaxArgMax<TIn, TOut>(Func<TIn, TOut> func, IEnumerable<TIn> inputs) where TOut : IComparable
+        {
+            var it = inputs.GetEnumerator();
+            it.MoveNext();
+            var maxArg = it.Current;
+            var maxValue = func(it.Current);
+            while (it.MoveNext())
+            {
+                var currentArg = it.Current;
+                var currentValue = func(it.Current);
+                if (maxValue.CompareTo(currentValue) < 0)
+                {
+                    maxArg = currentArg;
+                    maxValue = currentValue;
+                }
+            }
+
+            return Tuple.Create(maxArg, maxValue);
         }
     }