Free cookie consent management tool by TermsFeed Policy Generator

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

Last change on this file since 17618 was 17618, checked in by dleko, 4 years ago

#2825 Add storable generation field.
Add StorableType attribute to Solution class.

File size: 21.0 KB
Line 
1using System;
2using System.Collections.Generic;
3using System.Linq;
4using System.Threading;
5using HEAL.Attic;
6using HeuristicLab.Common;
7using HeuristicLab.Core;
8using HeuristicLab.Data;
9using HeuristicLab.Encodings.RealVectorEncoding;
10using HeuristicLab.Optimization;
11using HeuristicLab.Parameters;
12using HeuristicLab.Random;
13
14namespace HeuristicLab.Algorithms.NSGA3
15{
16    /// <summary>
17    /// The Reference Point Based Non-dominated Sorting Genetic Algorithm III was introduced in Deb
18    /// et al. 2013. An Evolutionary Many-Objective Optimization Algorithm Using Reference Point
19    /// Based Non-dominated Sorting Approach. IEEE Transactions on Evolutionary Computation, 18(4),
20    /// pp. 577-601.
21    /// </summary>
22    [Item("NSGA-III", "The Reference Point Based Non-dominated Sorting Genetic Algorithm III was introduced in Deb et al. 2013. An Evolutionary Many-Objective Optimization Algorithm Using Reference Point Based Non-dominated Sorting Approach. IEEE Transactions on Evolutionary Computation, 18(4), pp. 577-601.")]
23    [Creatable(Category = CreatableAttribute.Categories.PopulationBasedAlgorithms, Priority = 136)]
24    [StorableType("07C745F7-A8A3-4F99-8B2C-F97E639F9AC3")]
25    public class NSGA3 : BasicAlgorithm
26    {
27        private const double EPSILON = 10e-6; // a tiny number that is greater than 0
28
29        public override bool SupportsPause => false;
30
31        #region ProblemProperties
32
33        public override Type ProblemType
34        {
35            get { return typeof(MultiObjectiveBasicProblem<RealVectorEncoding>); }
36        }
37
38        public new MultiObjectiveBasicProblem<RealVectorEncoding> Problem
39        {
40            get { return (MultiObjectiveBasicProblem<RealVectorEncoding>)base.Problem; }
41            set { base.Problem = value; }
42        }
43
44        #endregion ProblemProperties
45
46        #region Storable fields
47
48        [Storable]
49        private IRandom random;
50
51        [Storable]
52        private int generation;
53
54        [Storable]
55        private List<Solution> solutions;
56
57        #endregion Storable fields
58
59        #region ParameterAndResultsNames
60
61        // Parameter Names
62
63        private const string SeedName = "Seed";
64        private const string SetSeedRandomlyName = "SetSeedRandomly";
65        private const string PopulationSizeName = "PopulationSize";
66        private const string CrossoverProbabilityName = "CrossOverProbability";
67        private const string MutationProbabilityName = "MutationProbability";
68        private const string MaximumGenerationsName = "MaximumGenerations";
69        private const string DominateOnEqualQualitiesName = "DominateOnEqualQualities";
70
71        // Results Names
72
73        private const string GeneratedReferencePointsResultName = "Generated Reference Points";
74        private const string CurrentFrontResultName = "Pareto Front"; // Do not touch this
75
76        #endregion ParameterAndResultsNames
77
78        #region ParameterProperties
79
80        private IFixedValueParameter<IntValue> SeedParameter
81        {
82            get { return (IFixedValueParameter<IntValue>)Parameters[SeedName]; }
83        }
84
85        private IFixedValueParameter<BoolValue> SetSeedRandomlyParameter
86        {
87            get { return (IFixedValueParameter<BoolValue>)Parameters[SetSeedRandomlyName]; }
88        }
89
90        private IFixedValueParameter<IntValue> PopulationSizeParameter
91        {
92            get { return (IFixedValueParameter<IntValue>)Parameters[PopulationSizeName]; }
93        }
94
95        private IFixedValueParameter<PercentValue> CrossoverProbabilityParameter
96        {
97            get { return (IFixedValueParameter<PercentValue>)Parameters[CrossoverProbabilityName]; }
98        }
99
100        private IFixedValueParameter<PercentValue> MutationProbabilityParameter
101        {
102            get { return (IFixedValueParameter<PercentValue>)Parameters[MutationProbabilityName]; }
103        }
104
105        private IFixedValueParameter<IntValue> MaximumGenerationsParameter
106        {
107            get { return (IFixedValueParameter<IntValue>)Parameters[MaximumGenerationsName]; }
108        }
109
110        private IFixedValueParameter<BoolValue> DominateOnEqualQualitiesParameter
111        {
112            get { return (IFixedValueParameter<BoolValue>)Parameters[DominateOnEqualQualitiesName]; }
113        }
114
115        #endregion ParameterProperties
116
117        #region Properties
118
119        public IntValue Seed => SeedParameter.Value;
120
121        public BoolValue SetSeedRandomly => SetSeedRandomlyParameter.Value;
122
123        public IntValue PopulationSize => PopulationSizeParameter.Value;
124
125        public PercentValue CrossoverProbability => CrossoverProbabilityParameter.Value;
126
127        public PercentValue MutationProbability => MutationProbabilityParameter.Value;
128
129        public IntValue MaximumGenerations => MaximumGenerationsParameter.Value;
130
131        public BoolValue DominateOnEqualQualities => DominateOnEqualQualitiesParameter.Value;
132
133        public List<List<Solution>> Fronts { get; private set; }
134
135        public List<ReferencePoint> ReferencePoints { get; private set; }
136
137        #endregion Properties
138
139        #region ResultsProperties
140
141        public DoubleMatrix ResultsGeneratedReferencePoints
142        {
143            get { return (DoubleMatrix)Results[GeneratedReferencePointsResultName].Value; }
144            set { Results[GeneratedReferencePointsResultName].Value = value; }
145        }
146
147        public DoubleMatrix ResultsSolutions
148        {
149            get { return (DoubleMatrix)Results[CurrentFrontResultName].Value; }
150            set { Results[CurrentFrontResultName].Value = value; }
151        }
152
153        #endregion ResultsProperties
154
155        #region Constructors
156
157        public NSGA3() : base()
158        {
159            Parameters.Add(new FixedValueParameter<IntValue>(SeedName, "The random seed used to initialize the new pseudo random number generator.", new IntValue(0)));
160            Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyName, "True if the random seed should be set to a random value, otherwise false.", new BoolValue(true)));
161            Parameters.Add(new FixedValueParameter<IntValue>(PopulationSizeName, "The size of the population of Individuals.", new IntValue(100)));
162            Parameters.Add(new FixedValueParameter<PercentValue>(CrossoverProbabilityName, "The probability that the crossover operator is applied on two parents.", new PercentValue(0.9)));
163            Parameters.Add(new FixedValueParameter<PercentValue>(MutationProbabilityName, "The probability that the mutation operator is applied on a Individual.", new PercentValue(0.05)));
164            Parameters.Add(new FixedValueParameter<IntValue>(MaximumGenerationsName, "The maximum number of generations which should be processed.", new IntValue(1000)));
165            Parameters.Add(new FixedValueParameter<BoolValue>(DominateOnEqualQualitiesName, "Flag which determines wether Individuals with equal quality values should be treated as dominated.", new BoolValue(false)));
166        }
167
168        // Persistence uses this ctor to improve deserialization efficiency. If we would use the
169        // default ctor instead this would completely initialize the object (e.g. creating
170        // parameters) even though the data is later overwritten by the stored data.
171        [StorableConstructor]
172        public NSGA3(StorableConstructorFlag _) : base(_) { }
173
174        // Each clonable item must have a cloning ctor (deep cloning, the cloner is used to handle
175        // cyclic object references). Don't forget to call the cloning ctor of the base class
176        public NSGA3(NSGA3 original, Cloner cloner) : base(original, cloner)
177        {
178            // todo: don't forget to clone storable fields
179            random = cloner.Clone(original.random);
180            solutions = new List<Solution>(original.solutions?.Select(cloner.Clone));
181        }
182
183        public override IDeepCloneable Clone(Cloner cloner)
184        {
185            return new NSGA3(this, cloner);
186        }
187
188        #endregion Constructors
189
190        #region Initialization
191
192        protected override void Initialize(CancellationToken cancellationToken)
193        {
194            base.Initialize(cancellationToken);
195
196            InitFields();
197            InitReferencePoints();
198            InitResults();
199            Analyze();
200        }
201
202        private void InitFields()
203        {
204            random = new MersenneTwister();
205            generation = 0;
206            InitSolutions();
207        }
208
209        private void InitSolutions()
210        {
211            int minBound = 0;
212            int maxBound = 1;
213
214            // Initialise solutions
215            solutions = new List<Solution>(PopulationSize.Value);
216            for (int i = 0; i < PopulationSize.Value; i++)
217            {
218                RealVector randomRealVector = new RealVector(Problem.Encoding.Length, random, minBound, maxBound);
219
220                solutions.Add(new Solution(StorableConstructorFlag.Default)
221                {
222                    Chromosome = randomRealVector
223                });
224                solutions[i].Fitness = Evaluate(solutions[i].Chromosome);
225            }
226        }
227
228        private void InitReferencePoints()
229        {
230            // Generate reference points and add them to results
231            int nDiv = 5; // todo: figure out the correct number of divisions
232            ReferencePoints = ReferencePoint.GenerateReferencePoints(Problem.Encoding.Length, nDiv);
233            ResultsGeneratedReferencePoints = Utility.ConvertToDoubleMatrix(ReferencePoints);
234        }
235
236        private void InitResults()
237        {
238            Results.Add(new Result(GeneratedReferencePointsResultName, "The initially generated reference points", Utility.ConvertToDoubleMatrix(ReferencePoints)));
239            Results.Add(new Result(CurrentFrontResultName, "The Pareto Front", new DoubleMatrix()));
240        }
241
242        #endregion Initialization
243
244        #region Overriden Methods
245
246        protected override void Run(CancellationToken cancellationToken)
247        {
248            ReferencePoints = new List<ReferencePoint>(); // todo: use existing list of reference points
249            while (generation != MaximumGenerations.Value)
250            {
251                ToNextGeneration();
252                generation++;
253            }
254        }
255
256        #endregion Overriden Methods
257
258        #region Private Methods
259
260        private void Analyze()
261        {
262            ResultsSolutions = solutions.Select(s => s.Chromosome.ToArray()).ToMatrix();
263            Problem.Analyze(
264                solutions.Select(s => (Individual)new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, s.Chromosome) } })).ToArray(),
265                solutions.Select(s => s.Fitness).ToArray(),
266                Results,
267                random
268                );
269        }
270
271        /// <summary>
272        /// Returns the fitness of the given <paramref name="chromosome" /> by applying the Evaluate
273        /// method of the Problem.
274        /// </summary>
275        /// <param name="chromosome"></param>
276        /// <returns></returns>
277        private double[] Evaluate(RealVector chromosome)
278        {
279            return Problem.Evaluate(new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, chromosome) } }), random);
280        }
281
282        private List<Solution> ToNextGeneration()
283        {
284            List<Solution> st = new List<Solution>();
285            List<Solution> qt = Mutate(Recombine(solutions));
286            List<Solution> rt = Utility.Concat(solutions, qt);
287            List<Solution> nextGeneration;
288
289            // Do non-dominated sort
290            var qualities = Utility.ToFitnessMatrix(solutions);
291            // compute the pareto fronts using the DominationCalculator and discard the qualities
292            // part in the inner tuples
293            Fronts = DominationCalculator<Solution>.CalculateAllParetoFronts(rt.ToArray(), qualities, Problem.Maximization, out int[] rank, true)
294                .Select(list => new List<Solution>(list.Select(pair => pair.Item1))).ToList();
295
296            int i = 0;
297            List<Solution> lf = null; // last front to be included
298            while (i < Fronts.Count && st.Count < PopulationSize.Value)
299            {
300                lf = Fronts[i];
301                st = Utility.Concat(st, lf);
302                i++;
303            }
304
305            if (st.Count == PopulationSize.Value) // no selection needs to be done
306                nextGeneration = st;
307            else
308            {
309                int l = i - 1;
310                nextGeneration = new List<Solution>();
311                for (int f = 0; f < l; f++)
312                    nextGeneration = Utility.Concat(nextGeneration, Fronts[f]);
313                Normalize(st);
314                Associate();
315                throw new NotImplementedException();
316            }
317
318            throw new NotImplementedException();
319        }
320
321        private void Normalize(List<Solution> population)
322        {
323            // Find the ideal point
324            double[] idealPoint = new double[Problem.Encoding.Length];
325            for (int j = 0; j < Problem.Encoding.Length; j++)
326            {
327                // Compute ideal point
328                idealPoint[j] = Utility.Min(s => s.Fitness[j], solutions);
329
330                // Translate objectives
331                foreach (var solution in solutions)
332                    solution.Fitness[j] -= idealPoint[j];
333            }
334
335            // Find the extreme points
336            Solution[] extremePoints = new Solution[Problem.Encoding.Length];
337            for (int j = 0; j < Problem.Encoding.Length; j++)
338            {
339                // Compute extreme points
340                double[] weights = new double[Problem.Encoding.Length];
341                for (int i = 0; i < Problem.Encoding.Length; i++) weights[i] = EPSILON;
342                weights[j] = 1;
343                double func(Solution s) => ASF(s.Fitness, weights);
344                extremePoints[j] = Utility.ArgMin(func, solutions);
345            }
346
347            // Compute intercepts
348            List<double> intercepts = GetIntercepts(extremePoints.ToList());
349
350            // Normalize objectives
351            NormalizeObjectives(intercepts, idealPoint);
352
353            // Associate reference points to solutions
354            Associate();
355        }
356
357        private void NormalizeObjectives(List<double> intercepts, double[] idealPoint)
358        {
359            for (int f = 0; f < Fronts.Count; f++)
360            {
361                foreach (var solution in Fronts[f])
362                {
363                    for (int i = 0; i < Problem.Encoding.Length; i++)
364                    {
365                        if (Math.Abs(intercepts[i] - idealPoint[i]) > EPSILON)
366                        {
367                            solution.Fitness[i] = solution.Fitness[i] / (intercepts[i] - idealPoint[i]);
368                        }
369                        else
370                        {
371                            solution.Fitness[i] = solution.Fitness[i] / EPSILON;
372                        }
373                    }
374                }
375            }
376        }
377
378        private void Associate()
379        {
380            for (int f = 0; f < Fronts.Count; f++)
381            {
382                foreach (var solution in Fronts[f])
383                {
384                    // find reference point for which the perpendicular distance to the current
385                    // solution is the lowest
386                    var rpAndDist = Utility.MinArgMin(rp => GetPerpendicularDistance(rp.Values, solution.Fitness), ReferencePoints);
387
388                    //// todo: use ArgMin here
389                    //int min_rp = -1;
390                    //double min_dist = double.MaxValue;
391                    //for (int r = 0; r < referencePoints.Count; r++)
392                    //{
393                    //    double d = GetPerpendicularDistance(referencePoints[r].Values, solution.Fitness);
394                    //    if (d < min_dist)
395                    //    {
396                    //        min_dist = d;
397                    //        min_rp = r;
398                    //    }
399                    //}
400                    if (f + 1 != Fronts.Count)
401                    {
402                        // Todo: Add member for reference point on index min_rp
403                        throw new NotImplementedException();
404                    }
405                    else
406                    {
407                        // Todo: Add potential member for reference point on index min_rp
408                        throw new NotImplementedException();
409                    }
410                }
411            }
412        }
413
414        private double GetPerpendicularDistance(double[] values, double[] fitness)
415        {
416            double numerator = 0;
417            double denominator = 0;
418            for (int i = 0; i < values.Length; i++)
419            {
420                numerator += values[i] * fitness[i];
421                denominator += Math.Pow(values[i], 2);
422            }
423            double k = numerator / denominator;
424
425            double d = 0;
426            for (int i = 0; i < values.Length; i++)
427            {
428                d += Math.Pow(k * values[i] - fitness[i], 2);
429            }
430            return Math.Sqrt(d);
431        }
432
433        private double ASF(double[] x, double[] weight)
434        {
435            List<int> dimensions = new List<int>();
436            for (int i = 0; i < Problem.Encoding.Length; i++) dimensions.Add(i);
437            double f(int dim) => x[dim] / weight[dim];
438            return Utility.Max(f, dimensions);
439        }
440
441        private List<double> GetIntercepts(List<Solution> extremePoints)
442        {
443            // Check whether there are duplicate extreme points. This might happen but the original
444            // paper does not mention how to deal with it.
445            bool duplicate = false;
446            for (int i = 0; !duplicate && i < extremePoints.Count; i++)
447            {
448                for (int j = i + 1; !duplicate && j < extremePoints.Count; j++)
449                {
450                    // maybe todo: override Equals method of solution?
451                    duplicate = extremePoints[i].Equals(extremePoints[j]);
452                }
453            }
454
455            List<double> intercepts = new List<double>();
456
457            if (duplicate)
458            { // cannot construct the unique hyperplane (this is a casual method to deal with the condition)
459                for (int f = 0; f < Problem.Encoding.Length; f++)
460                {
461                    // extreme_points[f] stands for the individual with the largest value of
462                    // objective f
463                    intercepts.Add(extremePoints[f].Fitness[f]);
464                }
465            }
466            else
467            {
468                // Find the equation of the hyperplane
469                List<double> b = new List<double>(); //(pop[0].objs().size(), 1.0);
470                for (int i = 0; i < Problem.Encoding.Length; i++)
471                {
472                    b.Add(1.0);
473                }
474
475                List<List<double>> a = new List<List<double>>();
476                foreach (Solution s in extremePoints)
477                {
478                    List<double> aux = new List<double>();
479                    for (int i = 0; i < Problem.Encoding.Length; i++)
480                        aux.Add(s.Fitness[i]);
481                    a.Add(aux);
482                }
483                List<double> x = GaussianElimination(a, b);
484
485                // Find intercepts
486                for (int f = 0; f < Problem.Encoding.Length; f++)
487                {
488                    intercepts.Add(1.0 / x[f]);
489                }
490            }
491
492            return intercepts;
493        }
494
495        private List<double> GaussianElimination(List<List<double>> a, List<double> b)
496        {
497            List<double> x = new List<double>();
498
499            int n = a.Count;
500            for (int i = 0; i < n; i++)
501                a[i].Add(b[i]);
502
503            for (int @base = 0; @base < n - 1; @base++)
504                for (int target = @base + 1; target < n; target++)
505                {
506                    double ratio = a[target][@base] / a[@base][@base];
507                    for (int term = 0; term < a[@base].Count; term++)
508                        a[target][term] = a[target][term] - a[@base][term] * ratio;
509                }
510
511            for (int i = 0; i < n; i++)
512                x.Add(0.0);
513
514            for (int i = n - 1; i >= 0; i--)
515            {
516                for (int known = i + 1; known < n; known++)
517                    a[i][n] = a[i][n] - a[i][known] * x[known];
518                x[i] = a[i][n] / a[i][i];
519            }
520
521            return x;
522        }
523
524        private List<Solution> Recombine(List<Solution> solutions)
525        {
526            throw new NotImplementedException();
527        }
528
529        private List<Solution> Mutate(List<Solution> solutions)
530        {
531            throw new NotImplementedException();
532        }
533
534        #endregion Private Methods
535    }
536}
Note: See TracBrowser for help on using the repository browser.