source: branches/MOCMAEvolutionStrategy/HeuristicLab.Algorithms.MOCMAEvolutionStrategy/3.3/MOCMAEvolutionStrategy.cs @ 15045

Last change on this file since 15045 was 15045, checked in by bwerth, 4 years ago

#2592 improvements/changes as requested in review comment

File size: 30.1 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
4 * and the BEACON Center for the Study of Evolution in Action.
5 *
6 * This file is part of HeuristicLab.
7 *
8 * HeuristicLab is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * HeuristicLab is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
20 */
21#endregion
22
23using System;
24using System.Collections.Generic;
25using System.Linq;
26using System.Threading;
27using HeuristicLab.Analysis;
28using HeuristicLab.Common;
29using HeuristicLab.Core;
30using HeuristicLab.Data;
31using HeuristicLab.Encodings.RealVectorEncoding;
32using HeuristicLab.Optimization;
33using HeuristicLab.Parameters;
34using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
35using HeuristicLab.Problems.TestFunctions.MultiObjective;
36using HeuristicLab.Random;
37
38namespace HeuristicLab.Algorithms.MOCMAEvolutionStrategy {
39  [Item("MOCMA Evolution Strategy (MOCMAES)", "A multi objective evolution strategy based on covariance matrix adaptation. Code is based on 'Covariance Matrix Adaptation for Multi - objective Optimization' by Igel, Hansen and Roth")]
40  [Creatable(CreatableAttribute.Categories.PopulationBasedAlgorithms, Priority = 210)]
41  [StorableClass]
42  [System.Runtime.InteropServices.Guid("5AC20A69-BBBF-4153-B57D-3EAF92DC505E")]
43  public class MOCMAEvolutionStrategy : BasicAlgorithm {
44    public override Type ProblemType
45    {
46      get { return typeof(MultiObjectiveBasicProblem<RealVectorEncoding>); }
47    }
48    public new MultiObjectiveBasicProblem<RealVectorEncoding> Problem
49    {
50      get { return (MultiObjectiveBasicProblem<RealVectorEncoding>)base.Problem; }
51      set { base.Problem = value; }
52    }
53    public override bool SupportsPause
54    {
55      get { return true; }
56    }
57
58    #region storable fields
59    [Storable]
60    private IRandom random = new MersenneTwister();
61    [Storable]
62    private NormalDistributedRandom gauss;
63    [Storable]
64    private Individual[] solutions;
65    [Storable]
66    private double stepSizeLearningRate; //=cp learning rate in [0,1]
67    [Storable]
68    private double stepSizeDampeningFactor; //d
69    [Storable]
70    private double targetSuccessProbability;// p^target_succ
71    [Storable]
72    private double evolutionPathLearningRate;//cc
73    [Storable]
74    private double covarianceMatrixLearningRate;//ccov
75    [Storable]
76    private double covarianceMatrixUnlearningRate;
77    [Storable]
78    private double successThreshold; //ptresh
79    #endregion
80
81    #region ParameterNames
82    private const string MaximumRuntimeName = "Maximum Runtime";
83    private const string SeedName = "Seed";
84    private const string SetSeedRandomlyName = "SetSeedRandomly";
85    private const string PopulationSizeName = "PopulationSize";
86    private const string MaximumGenerationsName = "MaximumGenerations";
87    private const string MaximumEvaluatedSolutionsName = "MaximumEvaluatedSolutions";
88    private const string InitialSigmaName = "InitialSigma";
89    private const string IndicatorName = "Indicator";
90
91    private const string EvaluationsResultName = "Evaluations";
92    private const string IterationsResultName = "Generations";
93    private const string TimetableResultName = "Timetable";
94    private const string HypervolumeResultName = "Hypervolume";
95    private const string GenerationalDistanceResultName = "Generational Distance";
96    private const string InvertedGenerationalDistanceResultName = "Inverted Generational Distance";
97    private const string CrowdingResultName = "Crowding";
98    private const string SpacingResultName = "Spacing";
99    private const string CurrentFrontResultName = "Pareto Front";
100    private const string BestHypervolumeResultName = "Best Hypervolume";
101    private const string BestKnownHypervolumeResultName = "Best known hypervolume";
102    private const string DifferenceToBestKnownHypervolumeResultName = "Absolute Distance to BestKnownHypervolume";
103    private const string ScatterPlotResultName = "ScatterPlot";
104    #endregion
105
106    #region ParameterProperties
107    public IFixedValueParameter<IntValue> MaximumRuntimeParameter
108    {
109      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumRuntimeName]; }
110    }
111    public IFixedValueParameter<IntValue> SeedParameter
112    {
113      get { return (IFixedValueParameter<IntValue>)Parameters[SeedName]; }
114    }
115    public FixedValueParameter<BoolValue> SetSeedRandomlyParameter
116    {
117      get { return (FixedValueParameter<BoolValue>)Parameters[SetSeedRandomlyName]; }
118    }
119    public IFixedValueParameter<IntValue> PopulationSizeParameter
120    {
121      get { return (IFixedValueParameter<IntValue>)Parameters[PopulationSizeName]; }
122    }
123    public IFixedValueParameter<IntValue> MaximumGenerationsParameter
124    {
125      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumGenerationsName]; }
126    }
127    public IFixedValueParameter<IntValue> MaximumEvaluatedSolutionsParameter
128    {
129      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumEvaluatedSolutionsName]; }
130    }
131    public IValueParameter<DoubleArray> InitialSigmaParameter
132    {
133      get { return (IValueParameter<DoubleArray>)Parameters[InitialSigmaName]; }
134    }
135    public IConstrainedValueParameter<IIndicator> IndicatorParameter
136    {
137      get { return (IConstrainedValueParameter<IIndicator>)Parameters[IndicatorName]; }
138    }
139    #endregion
140
141    #region Properties
142    public int MaximumRuntime
143    {
144      get { return MaximumRuntimeParameter.Value.Value; }
145      set { MaximumRuntimeParameter.Value.Value = value; }
146    }
147    public int Seed
148    {
149      get { return SeedParameter.Value.Value; }
150      set { SeedParameter.Value.Value = value; }
151    }
152    public bool SetSeedRandomly
153    {
154      get { return SetSeedRandomlyParameter.Value.Value; }
155      set { SetSeedRandomlyParameter.Value.Value = value; }
156    }
157    public int PopulationSize
158    {
159      get { return PopulationSizeParameter.Value.Value; }
160      set { PopulationSizeParameter.Value.Value = value; }
161    }
162    public int MaximumGenerations
163    {
164      get { return MaximumGenerationsParameter.Value.Value; }
165      set { MaximumGenerationsParameter.Value.Value = value; }
166    }
167    public int MaximumEvaluatedSolutions
168    {
169      get { return MaximumEvaluatedSolutionsParameter.Value.Value; }
170      set { MaximumEvaluatedSolutionsParameter.Value.Value = value; }
171    }
172    public DoubleArray InitialSigma
173    {
174      get { return InitialSigmaParameter.Value; }
175      set { InitialSigmaParameter.Value = value; }
176    }
177    public IIndicator Indicator
178    {
179      get { return IndicatorParameter.Value; }
180      set { IndicatorParameter.Value = value; }
181    }
182
183    public double StepSizeLearningRate { get { return stepSizeLearningRate; } }
184    public double StepSizeDampeningFactor { get { return stepSizeDampeningFactor; } }
185    public double TargetSuccessProbability { get { return targetSuccessProbability; } }
186    public double EvolutionPathLearningRate { get { return evolutionPathLearningRate; } }
187    public double CovarianceMatrixLearningRate { get { return covarianceMatrixLearningRate; } }
188    public double CovarianceMatrixUnlearningRate { get { return covarianceMatrixUnlearningRate; } }
189    public double SuccessThreshold { get { return successThreshold; } }
190    #endregion
191
192    #region ResultsProperties
193    private int ResultsEvaluations
194    {
195      get { return ((IntValue)Results[EvaluationsResultName].Value).Value; }
196      set { ((IntValue)Results[EvaluationsResultName].Value).Value = value; }
197    }
198    private int ResultsIterations
199    {
200      get { return ((IntValue)Results[IterationsResultName].Value).Value; }
201      set { ((IntValue)Results[IterationsResultName].Value).Value = value; }
202    }
203    #region Datatable
204    private DataTable ResultsQualities
205    {
206      get { return (DataTable)Results[TimetableResultName].Value; }
207    }
208    private DataRow ResultsBestHypervolumeDataLine
209    {
210      get { return ResultsQualities.Rows[BestHypervolumeResultName]; }
211    }
212    private DataRow ResultsHypervolumeDataLine
213    {
214      get { return ResultsQualities.Rows[HypervolumeResultName]; }
215    }
216    private DataRow ResultsGenerationalDistanceDataLine
217    {
218      get { return ResultsQualities.Rows[GenerationalDistanceResultName]; }
219    }
220    private DataRow ResultsInvertedGenerationalDistanceDataLine
221    {
222      get { return ResultsQualities.Rows[InvertedGenerationalDistanceResultName]; }
223    }
224    private DataRow ResultsCrowdingDataLine
225    {
226      get { return ResultsQualities.Rows[CrowdingResultName]; }
227    }
228    private DataRow ResultsSpacingDataLine
229    {
230      get { return ResultsQualities.Rows[SpacingResultName]; }
231    }
232    private DataRow ResultsHypervolumeDifferenceDataLine
233    {
234      get { return ResultsQualities.Rows[DifferenceToBestKnownHypervolumeResultName]; }
235    }
236    #endregion
237    //QualityIndicators
238    private double ResultsHypervolume
239    {
240      get { return ((DoubleValue)Results[HypervolumeResultName].Value).Value; }
241      set { ((DoubleValue)Results[HypervolumeResultName].Value).Value = value; }
242    }
243    private double ResultsGenerationalDistance
244    {
245      get { return ((DoubleValue)Results[GenerationalDistanceResultName].Value).Value; }
246      set { ((DoubleValue)Results[GenerationalDistanceResultName].Value).Value = value; }
247    }
248    private double ResultsInvertedGenerationalDistance
249    {
250      get { return ((DoubleValue)Results[InvertedGenerationalDistanceResultName].Value).Value; }
251      set { ((DoubleValue)Results[InvertedGenerationalDistanceResultName].Value).Value = value; }
252    }
253    private double ResultsCrowding
254    {
255      get { return ((DoubleValue)Results[CrowdingResultName].Value).Value; }
256      set { ((DoubleValue)Results[CrowdingResultName].Value).Value = value; }
257    }
258    private double ResultsSpacing
259    {
260      get { return ((DoubleValue)Results[SpacingResultName].Value).Value; }
261      set { ((DoubleValue)Results[SpacingResultName].Value).Value = value; }
262    }
263    private double ResultsBestHypervolume
264    {
265      get { return ((DoubleValue)Results[BestHypervolumeResultName].Value).Value; }
266      set { ((DoubleValue)Results[BestHypervolumeResultName].Value).Value = value; }
267    }
268    private double ResultsBestKnownHypervolume
269    {
270      get { return ((DoubleValue)Results[BestKnownHypervolumeResultName].Value).Value; }
271      set { ((DoubleValue)Results[BestKnownHypervolumeResultName].Value).Value = value; }
272    }
273    private double ResultsDifferenceBestKnownHypervolume
274    {
275      get { return ((DoubleValue)Results[DifferenceToBestKnownHypervolumeResultName].Value).Value; }
276      set { ((DoubleValue)Results[DifferenceToBestKnownHypervolumeResultName].Value).Value = value; }
277
278    }
279    //Solutions
280    private DoubleMatrix ResultsSolutions
281    {
282      get { return (DoubleMatrix)Results[CurrentFrontResultName].Value; }
283      set { Results[CurrentFrontResultName].Value = value; }
284    }
285    private ScatterPlotContent ResultsScatterPlot
286    {
287      get { return (ScatterPlotContent)Results[ScatterPlotResultName].Value; }
288      set { Results[ScatterPlotResultName].Value = value; }
289    }
290    #endregion
291
292    #region Constructors
293    public MOCMAEvolutionStrategy() {
294      Parameters.Add(new FixedValueParameter<IntValue>(MaximumRuntimeName, "The maximum runtime in seconds after which the algorithm stops. Use -1 to specify no limit for the runtime", new IntValue(3600)));
295      Parameters.Add(new FixedValueParameter<IntValue>(SeedName, "The random seed used to initialize the new pseudo random number generator.", new IntValue(0)));
296      Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyName, "True if the random seed should be set to a random value, otherwise false.", new BoolValue(true)));
297      Parameters.Add(new FixedValueParameter<IntValue>(PopulationSizeName, "λ (lambda) - the size of the offspring population.", new IntValue(20)));
298      Parameters.Add(new ValueParameter<DoubleArray>(InitialSigmaName, "The initial sigma can be a single value or a value for each dimension. All values need to be > 0.", new DoubleArray(new[] { 0.5 })));
299      Parameters.Add(new FixedValueParameter<IntValue>(MaximumGenerationsName, "The maximum number of generations which should be processed.", new IntValue(1000)));
300      Parameters.Add(new FixedValueParameter<IntValue>(MaximumEvaluatedSolutionsName, "The maximum number of evaluated solutions that should be computed.", new IntValue(int.MaxValue)));
301      var set = new ItemSet<IIndicator> { new HypervolumeIndicator(), new CrowdingIndicator(), new MinimalDistanceIndicator() };
302      Parameters.Add(new ConstrainedValueParameter<IIndicator>(IndicatorName, "The selection mechanism on non-dominated solutions", set, set.First()));
303    }
304
305    [StorableConstructor]
306    protected MOCMAEvolutionStrategy(bool deserializing) : base(deserializing) { }
307
308    protected MOCMAEvolutionStrategy(MOCMAEvolutionStrategy original, Cloner cloner) : base(original, cloner) {
309      random = cloner.Clone(original.random);
310      gauss = cloner.Clone(original.gauss);
311      solutions = original.solutions.Select(cloner.Clone).ToArray();
312      stepSizeLearningRate = original.stepSizeLearningRate;
313      stepSizeDampeningFactor = original.stepSizeDampeningFactor;
314      targetSuccessProbability = original.targetSuccessProbability;
315      evolutionPathLearningRate = original.evolutionPathLearningRate;
316      covarianceMatrixLearningRate = original.covarianceMatrixLearningRate;
317      covarianceMatrixUnlearningRate = original.covarianceMatrixUnlearningRate;
318      successThreshold = original.successThreshold;
319    }
320
321    public override IDeepCloneable Clone(Cloner cloner) { return new MOCMAEvolutionStrategy(this, cloner); }
322    #endregion
323
324    #region Initialization
325    protected override void Initialize(CancellationToken cancellationToken) {
326      if (SetSeedRandomly) Seed = new System.Random().Next();
327      random.Reset(Seed);
328      gauss = new NormalDistributedRandom(random, 0, 1);
329
330      InitResults();
331      InitStrategy();
332      InitSolutions();
333      Analyze();
334
335      ResultsIterations = 1;
336      cancellationToken.ThrowIfCancellationRequested();
337    }
338    private Individual InitializeIndividual(RealVector x) {
339      var zeros = new RealVector(x.Length);
340      var c = new double[x.Length, x.Length];
341      var sigma = InitialSigma.Max();
342      for (var i = 0; i < x.Length; i++) {
343        var d = InitialSigma[i % InitialSigma.Length] / sigma;
344        c[i, i] = d * d;
345      }
346      return new Individual(x, targetSuccessProbability, sigma, zeros, c, this);
347    }
348    private void InitSolutions() {
349      solutions = new Individual[PopulationSize];
350      for (var i = 0; i < PopulationSize; i++) {
351        var x = new RealVector(Problem.Encoding.Length); // Uniform distibution in all dimensions assumed.
352        var bounds = Problem.Encoding.Bounds;
353        for (var j = 0; j < Problem.Encoding.Length; j++) {
354          var dim = j % bounds.Rows;
355          x[j] = random.NextDouble() * (bounds[dim, 1] - bounds[dim, 0]) + bounds[dim, 0];
356        }
357        solutions[i] = InitializeIndividual(x);
358        PenalizeEvaluate(solutions[i]);
359      }
360    }
361    private void InitStrategy() {
362      const int lambda = 1;
363      double n = Problem.Encoding.Length;
364      targetSuccessProbability = 1.0 / (5.0 + Math.Sqrt(lambda) / 2.0);
365      stepSizeDampeningFactor = 1.0 + n / (2.0 * lambda);
366      stepSizeLearningRate = targetSuccessProbability * lambda / (2.0 + targetSuccessProbability * lambda);
367      evolutionPathLearningRate = 2.0 / (n + 2.0);
368      covarianceMatrixLearningRate = 2.0 / (n * n + 6.0);
369      covarianceMatrixUnlearningRate = 0.4 / (Math.Pow(n, 1.6) + 1);
370      successThreshold = 0.44;
371    }
372    private void InitResults() {
373      Results.Add(new Result(IterationsResultName, "The number of gererations evaluated", new IntValue(0)));
374      Results.Add(new Result(EvaluationsResultName, "The number of function evaltions performed", new IntValue(0)));
375      Results.Add(new Result(HypervolumeResultName, "The hypervolume of the current front considering the Referencepoint defined in the Problem", new DoubleValue(0.0)));
376      Results.Add(new Result(BestHypervolumeResultName, "The best hypervolume of the current run considering the Referencepoint defined in the Problem", new DoubleValue(0.0)));
377      Results.Add(new Result(BestKnownHypervolumeResultName, "The best knwon hypervolume considering the Referencepoint defined in the Problem", new DoubleValue(double.NaN)));
378      Results.Add(new Result(DifferenceToBestKnownHypervolumeResultName, "The difference between the current and the best known hypervolume", new DoubleValue(double.NaN)));
379      Results.Add(new Result(GenerationalDistanceResultName, "The generational distance to an optimal pareto front defined in the Problem", new DoubleValue(double.NaN)));
380      Results.Add(new Result(InvertedGenerationalDistanceResultName, "The inverted generational distance to an optimal pareto front defined in the Problem", new DoubleValue(double.NaN)));
381      Results.Add(new Result(CrowdingResultName, "The average crowding value for the current front (excluding infinities)", new DoubleValue(0.0)));
382      Results.Add(new Result(SpacingResultName, "The spacing for the current front (excluding infinities)", new DoubleValue(0.0)));
383
384      var table = new DataTable("QualityIndicators");
385      table.Rows.Add(new DataRow(BestHypervolumeResultName));
386      table.Rows.Add(new DataRow(HypervolumeResultName));
387      table.Rows.Add(new DataRow(CrowdingResultName));
388      table.Rows.Add(new DataRow(GenerationalDistanceResultName));
389      table.Rows.Add(new DataRow(InvertedGenerationalDistanceResultName));
390      table.Rows.Add(new DataRow(DifferenceToBestKnownHypervolumeResultName));
391      table.Rows.Add(new DataRow(SpacingResultName));
392      Results.Add(new Result(TimetableResultName, "Different quality meassures in a timeseries", table));
393      Results.Add(new Result(CurrentFrontResultName, "The current front", new DoubleMatrix()));
394      Results.Add(new Result(ScatterPlotResultName, "A scatterplot displaying the evaluated solutions and (if available) the analytically optimal front", new ScatterPlotContent(null, null, null, 2)));
395
396      var problem = Problem as MultiObjectiveTestFunctionProblem;
397      if (problem == null) return;
398      if (problem.BestKnownFront != null) {
399        ResultsBestKnownHypervolume = Hypervolume.Calculate(problem.BestKnownFront.ToJaggedArray(), problem.TestFunction.ReferencePoint(problem.Objectives), Problem.Maximization);
400        ResultsDifferenceBestKnownHypervolume = ResultsBestKnownHypervolume;
401      }
402      //TODO? move FrontScatterPlotContent partially? to MultiobjectiveTestProblem?
403      ResultsScatterPlot = new ScatterPlotContent(new double[0][], new double[0][], problem.BestKnownFront.ToJaggedArray(), problem.Objectives);
404    }
405    #endregion
406
407    #region Mainloop
408    protected override void Run(CancellationToken cancellationToken) {
409      while (ResultsIterations < MaximumGenerations) {
410        try {
411          Iterate();
412          ResultsIterations++;
413          cancellationToken.ThrowIfCancellationRequested();
414        }
415        finally {
416          Analyze();
417        }
418      }
419    }
420    private void Iterate() {
421      var offspring = solutions.Select(i => {
422        var o = new Individual(i);
423        o.Mutate(gauss);
424        PenalizeEvaluate(o);
425        return o;
426      });
427      var parents = solutions.Concat(offspring).ToArray();
428      SelectParents(parents, solutions.Length);
429      UpdatePopulation(parents);
430    }
431    protected override void OnExecutionTimeChanged() {
432      base.OnExecutionTimeChanged();
433      if (CancellationTokenSource == null) return;
434      if (MaximumRuntime == -1) return;
435      if (ExecutionTime.TotalSeconds > MaximumRuntime) CancellationTokenSource.Cancel();
436    }
437    #endregion
438
439    #region Evaluation
440    private void PenalizeEvaluate(Individual individual) {
441      if (IsFeasable(individual.Mean)) {
442        individual.Fitness = Evaluate(individual.Mean);
443        individual.PenalizedFitness = individual.Fitness;
444      } else {
445        var t = ClosestFeasible(individual.Mean);
446        individual.Fitness = Evaluate(t);
447        individual.PenalizedFitness = Penalize(individual.Mean, t, individual.Fitness);
448      }
449    }
450    private double[] Evaluate(RealVector x) {
451      var res = Problem.Evaluate(new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, x) } }), random);
452      ResultsEvaluations++;
453      return res;
454    }
455    private double[] Penalize(RealVector x, RealVector t, IEnumerable<double> fitness) {
456      var penalty = x.Zip(t, (a, b) => (a - b) * (a - b)).Sum() * 1E-6;
457      return fitness.Select((v, i) => Problem.Maximization[i] ? v - penalty : v + penalty).ToArray();
458    }
459    private RealVector ClosestFeasible(RealVector x) {
460      var bounds = Problem.Encoding.Bounds;
461      var r = new RealVector(x.Length);
462      for (var i = 0; i < x.Length; i++) {
463        var dim = i % bounds.Rows;
464        r[i] = Math.Min(Math.Max(bounds[dim, 0], x[i]), bounds[dim, 1]);
465      }
466      return r;
467    }
468    private bool IsFeasable(RealVector offspring) {
469      var bounds = Problem.Encoding.Bounds;
470      for (var i = 0; i < offspring.Length; i++) {
471        var dim = i % bounds.Rows;
472        if (bounds[dim, 0] > offspring[i] || offspring[i] > bounds[dim, 1]) return false;
473      }
474      return true;
475    }
476    #endregion
477
478    private void SelectParents(IReadOnlyList<Individual> parents, int length) {
479      //perform a nondominated sort to assign the rank to every element
480      var fronts = NonDominatedSort(parents);
481
482      //deselect the highest rank fronts until we would end up with less or equal mu elements
483      var rank = fronts.Count - 1;
484      var popSize = parents.Count;
485      while (popSize - fronts[rank].Count >= length) {
486        var front = fronts[rank];
487        foreach (var i in front) i.Selected = false;
488        popSize -= front.Count;
489        rank--;
490      }
491
492      //now use the indicator to deselect the approximatingly worst elements of the last selected front
493      var front1 = fronts[rank].OrderBy(x => x.PenalizedFitness[0]).ToList();
494      for (; popSize > length; popSize--) {
495        var lc = Indicator.LeastContributer(front1, Problem);
496        front1[lc].Selected = false;
497        front1.Swap(lc, front1.Count - 1);
498        front1.RemoveAt(front1.Count - 1);
499      }
500    }
501
502    private void UpdatePopulation(IReadOnlyList<Individual> parents) {
503      foreach (var p in parents.Skip(solutions.Length).Where(i => i.Selected))
504        p.UpdateAsOffspring();
505
506      for (var i = 0; i < solutions.Length; i++)
507        if (parents[i].Selected)
508          parents[i].UpdateAsParent(parents[i + solutions.Length].Selected);
509
510      solutions = parents.Where(p => p.Selected).ToArray();
511    }
512
513    private void Analyze() {
514      //TODO? move FrontScatterPlotContent partially to MultiobjectiveTestProblem
515      ResultsScatterPlot = new ScatterPlotContent(solutions.Select(x => x.Fitness).ToArray(), solutions.Select(x => x.Mean.ToArray()).ToArray(), ResultsScatterPlot.ParetoFront, ResultsScatterPlot.Objectives);
516
517      ResultsSolutions = solutions.Select(x => x.Mean.ToArray()).ToMatrix();
518
519      var problem = Problem as MultiObjectiveTestFunctionProblem;
520      if (problem == null) return;
521
522      var front = NonDominatedSelect.GetDominatingVectors(solutions.Select(x => x.Fitness), problem.ReferencePoint.CloneAsArray(), Problem.Maximization, true).ToArray();
523      if (front.Length == 0) return;
524      var bounds = problem.Bounds.CloneAsMatrix();
525      ResultsCrowding = Crowding.Calculate(front, bounds);
526      ResultsSpacing = Spacing.Calculate(front);
527      ResultsGenerationalDistance = problem.BestKnownFront != null ? GenerationalDistance.Calculate(front, problem.BestKnownFront.ToJaggedArray(), 1) : double.NaN;
528      ResultsInvertedGenerationalDistance = problem.BestKnownFront != null ? InvertedGenerationalDistance.Calculate(front, problem.BestKnownFront.ToJaggedArray(), 1) : double.NaN;
529      ResultsHypervolume = Hypervolume.Calculate(front, problem.ReferencePoint.CloneAsArray(), Problem.Maximization);
530      ResultsBestHypervolume = Math.Max(ResultsHypervolume, ResultsBestHypervolume);
531      ResultsDifferenceBestKnownHypervolume = ResultsBestKnownHypervolume - ResultsBestHypervolume;
532
533      ResultsBestHypervolumeDataLine.Values.Add(ResultsBestHypervolume);
534      ResultsHypervolumeDataLine.Values.Add(ResultsHypervolume);
535      ResultsCrowdingDataLine.Values.Add(ResultsCrowding);
536      ResultsGenerationalDistanceDataLine.Values.Add(ResultsGenerationalDistance);
537      ResultsInvertedGenerationalDistanceDataLine.Values.Add(ResultsInvertedGenerationalDistance);
538      ResultsSpacingDataLine.Values.Add(ResultsSpacing);
539      ResultsHypervolumeDifferenceDataLine.Values.Add(ResultsDifferenceBestKnownHypervolume);
540
541      Problem.Analyze(
542        solutions.Select(x => (Optimization.Individual)new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, x.Mean) } })).ToArray(),
543        solutions.Select(x => x.Fitness).ToArray(),
544        Results,
545        random);
546    }
547
548    #region FastNonDominatedSort
549    //blatantly stolen form HeuristicLab.Optimization.Operators.FastNonDominatedSort
550    //however: Operators.FastNonDominatedSort does not return ranked fronts => rerank after sorting would not be sensible
551
552    private enum DominationResult { Dominates, IsDominated, IsNonDominated };
553    private List<List<Individual>> NonDominatedSort(IReadOnlyList<Individual> individuals) {
554      const bool dominateOnEqualQualities = false;
555      var maximization = Problem.Maximization;
556      if (individuals == null) throw new InvalidOperationException(Name + ": No qualities found.");
557      var populationSize = individuals.Count;
558
559      var fronts = new List<List<Individual>>();
560      var dominatedScopes = new Dictionary<Individual, List<int>>();
561      var dominationCounter = new int[populationSize];
562
563      for (var pI = 0; pI < populationSize - 1; pI++) {
564        var p = individuals[pI];
565        List<int> dominatedScopesByp;
566        if (!dominatedScopes.TryGetValue(p, out dominatedScopesByp))
567          dominatedScopes[p] = dominatedScopesByp = new List<int>();
568        for (var qI = pI + 1; qI < populationSize; qI++) {
569          var test = Dominates(individuals[pI], individuals[qI], maximization, dominateOnEqualQualities);
570          if (test == DominationResult.Dominates) {
571            dominatedScopesByp.Add(qI);
572            dominationCounter[qI] += 1;
573          } else if (test == DominationResult.IsDominated) {
574            dominationCounter[pI] += 1;
575            if (!dominatedScopes.ContainsKey(individuals[qI]))
576              dominatedScopes.Add(individuals[qI], new List<int>());
577            dominatedScopes[individuals[qI]].Add(pI);
578          }
579          if (pI == populationSize - 2
580            && qI == populationSize - 1
581            && dominationCounter[qI] == 0) {
582            AddToFront(individuals[qI], fronts, 0);
583          }
584        }
585        if (dominationCounter[pI] == 0) {
586          AddToFront(p, fronts, 0);
587        }
588      }
589      var i = 0;
590      while (i < fronts.Count && fronts[i].Count > 0) {
591        var nextFront = new List<Individual>();
592        foreach (var p in fronts[i]) {
593          List<int> dominatedScopesByp;
594          if (!dominatedScopes.TryGetValue(p, out dominatedScopesByp)) continue;
595          foreach (var dominatedScope in dominatedScopesByp) {
596            dominationCounter[dominatedScope] -= 1;
597            if (dominationCounter[dominatedScope] != 0) continue;
598            nextFront.Add(individuals[dominatedScope]);
599          }
600        }
601        i += 1;
602        fronts.Add(nextFront);
603      }
604
605      for (i = 0; i < fronts.Count; i++) {
606        foreach (var p in fronts[i]) {
607          p.Rank = i;
608        }
609      }
610      return fronts;
611    }
612    private static void AddToFront(Individual p, IList<List<Individual>> fronts, int i) {
613      if (i == fronts.Count) fronts.Add(new List<Individual>());
614      fronts[i].Add(p);
615    }
616    private static DominationResult Dominates(Individual left, Individual right, bool[] maximizations, bool dominateOnEqualQualities) {
617      return Dominates(left.PenalizedFitness, right.PenalizedFitness, maximizations, dominateOnEqualQualities);
618    }
619    private static DominationResult Dominates(IReadOnlyList<double> left, IReadOnlyList<double> right, IReadOnlyList<bool> maximizations, bool dominateOnEqualQualities) {
620      //mkommend Caution: do not use LINQ.SequenceEqual for comparing the two quality arrays (left and right) due to performance reasons
621      if (dominateOnEqualQualities) {
622        var equal = true;
623        for (var i = 0; i < left.Count; i++) {
624          if (left[i] != right[i]) {
625            equal = false;
626            break;
627          }
628        }
629        if (equal) return DominationResult.Dominates;
630      }
631
632      bool leftIsBetter = false, rightIsBetter = false;
633      for (var i = 0; i < left.Count; i++) {
634        if (IsDominated(left[i], right[i], maximizations[i])) rightIsBetter = true;
635        else if (IsDominated(right[i], left[i], maximizations[i])) leftIsBetter = true;
636        if (leftIsBetter && rightIsBetter) break;
637      }
638
639      if (leftIsBetter && !rightIsBetter) return DominationResult.Dominates;
640      if (!leftIsBetter && rightIsBetter) return DominationResult.IsDominated;
641      return DominationResult.IsNonDominated;
642    }
643    private static bool IsDominated(double left, double right, bool maximization) {
644      return maximization && left < right
645        || !maximization && left > right;
646    }
647    #endregion
648  }
649}
Note: See TracBrowser for help on using the repository browser.