source: branches/MOCMAEvolutionStrategy/HeuristicLab.Algorithms.MOCMAEvolutionStrategy/3.3/MOCMASEvolutionStrategy.cs @ 14404

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2016 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 * and the BEACON Center for the Study of Evolution in Action.
 * 
 * This file is part of HeuristicLab.
 *
 * HeuristicLab is free software: you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation, either version 3 of the License, or
 * (at your option) any later version.
 *
 * HeuristicLab is distributed in the hope that it will be useful,
 * but WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 * GNU General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.
 */
#endregion

using System;
using System.Collections.Generic;
using System.Linq;
using System.Threading;
using HeuristicLab.Algorithms.CMAEvolutionStrategy;
using HeuristicLab.Analysis;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.RealVectorEncoding;
using HeuristicLab.Optimization;
using HeuristicLab.Parameters;
using HeuristicLab.Persistence.Default.CompositeSerializers.Storable;
using HeuristicLab.Problems.TestFunctions.MultiObjective;
using HeuristicLab.Random;

namespace HeuristicLab.Algorithms.MOCMAEvolutionStrategy {
  [Item("MOCMAS Evolution Strategy (MOCMASES)", "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")]
  [Creatable(CreatableAttribute.Categories.PopulationBasedAlgorithms, Priority = 210)]
  [StorableClass]
  public class MOCMASEvolutionStrategy : BasicAlgorithm {
    public override Type ProblemType
    {
      get { return typeof(MultiObjectiveTestFunctionProblem); }
    }
    public new MultiObjectiveTestFunctionProblem Problem
    {
      get { return (MultiObjectiveTestFunctionProblem)base.Problem; }
      set { base.Problem = value; }
    }
    #region internal variables
    private readonly IRandom random = new MersenneTwister();
    private NormalDistributedRandom gauss;
    private MOCMAESIndividual[] solutions;
    private MOCMAESParameters internals;
    #endregion

    #region ParameterNames
    private const string MaximumRuntimeName = "Maximum Runtime";
    private const string SeedName = "Seed";
    private const string SetSeedRandomlyName = "SetSeedRandomly";
    private const string PopulationSizeName = "PopulationSize";
    private const string MaximumGenerationsName = "MaximumGenerations";
    private const string MaximumEvaluatedSolutionsName = "MaximumEvaluatedSolutions";
    private const string InitialSigmaName = "InitialSigma";
    private const string IndicatorName = "Indicator";

    private const string EvaluationsResultName = "Evaluations";
    private const string IterationsResultName = "Generations";
    private const string TimetableResultName = "Timetable";
    private const string HypervolumeResultName = "Hypervolume";
    private const string GenerationalDistanceResultName = "Generational Distance";
    private const string InvertedGenerationalDistanceResultName = "Inverted Generational Distance";
    private const string CrowdingResultName = "Crowding";
    private const string SpacingResultName = "Spacing";
    private const string SolutionsResultName = "Pareto Front";
    private const string BestHypervolumeResultName = "Best Hypervolume";
    private const string BestKnownHypervolumeResultName = "Best known hypervolume";
    private const string DifferenceToBestKnownHypervolumeResultName = "Absolute Distance to BestKnownHypervolume";
    private const string ScatterPlotResultName = "ScatterPlot";

    #endregion

    #region ParameterProperties
    public IFixedValueParameter<IntValue> MaximumRuntimeParameter
    {
      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumRuntimeName]; }
    }
    public IFixedValueParameter<IntValue> SeedParameter
    {
      get { return (IFixedValueParameter<IntValue>)Parameters[SeedName]; }
    }
    public FixedValueParameter<BoolValue> SetSeedRandomlyParameter
    {
      get { return (FixedValueParameter<BoolValue>)Parameters[SetSeedRandomlyName]; }
    }
    private IFixedValueParameter<IntValue> PopulationSizeParameter
    {
      get { return (IFixedValueParameter<IntValue>)Parameters[PopulationSizeName]; }
    }
    private IFixedValueParameter<IntValue> MaximumGenerationsParameter
    {
      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumGenerationsName]; }
    }
    private IFixedValueParameter<IntValue> MaximumEvaluatedSolutionsParameter
    {
      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumEvaluatedSolutionsName]; }
    }
    private IFixedValueParameter<DoubleValue> InitialSigmaParameter
    {
      get { return (IFixedValueParameter<DoubleValue>)Parameters[InitialSigmaName]; }
    }
    private IConstrainedValueParameter<IIndicator> IndicatorParameter
    {
      get { return (IConstrainedValueParameter<IIndicator>)Parameters[IndicatorName]; }
    }
    #endregion

    #region Properties
    public int MaximumRuntime
    {
      get { return MaximumRuntimeParameter.Value.Value; }
      set { MaximumRuntimeParameter.Value.Value = value; }
    }
    public int Seed
    {
      get { return SeedParameter.Value.Value; }
      set { SeedParameter.Value.Value = value; }
    }
    public bool SetSeedRandomly
    {
      get { return SetSeedRandomlyParameter.Value.Value; }
      set { SetSeedRandomlyParameter.Value.Value = value; }
    }
    public int PopulationSize
    {
      get { return PopulationSizeParameter.Value.Value; }
      set { PopulationSizeParameter.Value.Value = value; }
    }
    public int MaximumGenerations
    {
      get { return MaximumGenerationsParameter.Value.Value; }
      set { MaximumGenerationsParameter.Value.Value = value; }
    }
    public int MaximumEvaluatedSolutions
    {
      get { return MaximumEvaluatedSolutionsParameter.Value.Value; }
      set { MaximumEvaluatedSolutionsParameter.Value.Value = value; }
    }
    public double InitialSigma
    {
      get { return InitialSigmaParameter.Value.Value; }
      set { InitialSigmaParameter.Value.Value = value; }
    }
    public IIndicator Indicator
    {
      get { return IndicatorParameter.Value; }
      set { IndicatorParameter.Value = value; }
    }


    #endregion

    #region ResultsProperties
    private int ResultsEvaluations
    {
      get { return ((IntValue)Results[EvaluationsResultName].Value).Value; }
      set { ((IntValue)Results[EvaluationsResultName].Value).Value = value; }
    }
    private int ResultsIterations
    {
      get { return ((IntValue)Results[IterationsResultName].Value).Value; }
      set { ((IntValue)Results[IterationsResultName].Value).Value = value; }
    }
    #region Datatable
    private DataTable ResultsQualities
    {
      get { return (DataTable)Results[TimetableResultName].Value; }
    }
    private DataRow ResultsBestHypervolumeDataLine
    {
      get { return ResultsQualities.Rows[BestHypervolumeResultName]; }
    }
    private DataRow ResultsHypervolumeDataLine
    {
      get { return ResultsQualities.Rows[HypervolumeResultName]; }
    }
    private DataRow ResultsGenerationalDistanceDataLine
    {
      get { return ResultsQualities.Rows[GenerationalDistanceResultName]; }
    }
    private DataRow ResultsInvertedGenerationalDistanceDataLine
    {
      get { return ResultsQualities.Rows[InvertedGenerationalDistanceResultName]; }
    }
    private DataRow ResultsCrowdingDataLine
    {
      get { return ResultsQualities.Rows[CrowdingResultName]; }
    }
    private DataRow ResultsSpacingDataLine
    {
      get { return ResultsQualities.Rows[SpacingResultName]; }
    }
    private DataRow ResultsHypervolumeDifferenceDataLine
    {
      get { return ResultsQualities.Rows[DifferenceToBestKnownHypervolumeResultName]; }
    }
    #endregion
    //QualityIndicators
    private double ResultsHypervolume
    {
      get { return ((DoubleValue)Results[HypervolumeResultName].Value).Value; }
      set { ((DoubleValue)Results[HypervolumeResultName].Value).Value = value; }
    }
    private double ResultsGenerationalDistance
    {
      get { return ((DoubleValue)Results[GenerationalDistanceResultName].Value).Value; }
      set { ((DoubleValue)Results[GenerationalDistanceResultName].Value).Value = value; }
    }
    private double ResultsInvertedGenerationalDistance
    {
      get { return ((DoubleValue)Results[InvertedGenerationalDistanceResultName].Value).Value; }
      set { ((DoubleValue)Results[InvertedGenerationalDistanceResultName].Value).Value = value; }
    }
    private double ResultsCrowding
    {
      get { return ((DoubleValue)Results[CrowdingResultName].Value).Value; }
      set { ((DoubleValue)Results[CrowdingResultName].Value).Value = value; }
    }
    private double ResultsSpacing
    {
      get { return ((DoubleValue)Results[SpacingResultName].Value).Value; }
      set { ((DoubleValue)Results[SpacingResultName].Value).Value = value; }
    }
    private double ResultsBestHypervolume
    {
      get { return ((DoubleValue)Results[BestHypervolumeResultName].Value).Value; }
      set { ((DoubleValue)Results[BestHypervolumeResultName].Value).Value = value; }
    }
    private double ResultsBestKnownHypervolume
    {
      get { return ((DoubleValue)Results[BestKnownHypervolumeResultName].Value).Value; }
      set { ((DoubleValue)Results[BestKnownHypervolumeResultName].Value).Value = value; }
    }
    private double ResultsDifferenceBestKnownHypervolume
    {
      get { return ((DoubleValue)Results[DifferenceToBestKnownHypervolumeResultName].Value).Value; }
      set { ((DoubleValue)Results[DifferenceToBestKnownHypervolumeResultName].Value).Value = value; }

    }
    //Solutions
    private DoubleMatrix ResultsSolutions
    {
      get { return ((DoubleMatrix)Results[SolutionsResultName].Value); }
      set { Results[SolutionsResultName].Value = value; }
    }
    private ScatterPlotContent ResultsScatterPlot
    {
      get { return ((ScatterPlotContent)Results[ScatterPlotResultName].Value); }
      set { Results[ScatterPlotResultName].Value = value; }
    }
    #endregion

    #region Constructors
    public MOCMASEvolutionStrategy() {
      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)));
      Parameters.Add(new FixedValueParameter<IntValue>(SeedName, "The random seed used to initialize the new pseudo random number generator.", new IntValue(0)));
      Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyName, "True if the random seed should be set to a random value, otherwise false.", new BoolValue(true)));
      Parameters.Add(new FixedValueParameter<IntValue>(PopulationSizeName, "λ (lambda) - the size of the offspring population.", new IntValue(20)));
      Parameters.Add(new FixedValueParameter<DoubleValue>(InitialSigmaName, "The initial sigma is a single value > 0.", new DoubleValue(0.5)));
      Parameters.Add(new FixedValueParameter<IntValue>(MaximumGenerationsName, "The maximum number of generations which should be processed.", new IntValue(1000)));
      Parameters.Add(new FixedValueParameter<IntValue>(MaximumEvaluatedSolutionsName, "The maximum number of evaluated solutions that should be computed.", new IntValue(int.MaxValue)));
      var set = new ItemSet<IIndicator> { new HypervolumeIndicator(), new CrowdingIndicator() };
      Parameters.Add(new ConstrainedValueParameter<IIndicator>(IndicatorName, "The selection mechanism on non-dominated solutions", set, set.First()));
    }

    [StorableConstructor]
    protected MOCMASEvolutionStrategy(bool deserializing) : base(deserializing) { }
    protected MOCMASEvolutionStrategy(MOCMASEvolutionStrategy original, Cloner cloner) : base(original, cloner) { }
    public override IDeepCloneable Clone(Cloner cloner) { return new MOCMASEvolutionStrategy(this, cloner); }
    #endregion

    #region Mainloop
    protected override void Run(CancellationToken cancellationToken) {
      // Set up the algorithm
      if (SetSeedRandomly) Seed = new System.Random().Next();
      random.Reset(Seed);
      gauss = new NormalDistributedRandom(random, 0, 1);

      InitResults();
      InitStrategy();
      InitSolutions();
      AnalyzeSolutions();

      // Loop until iteration limit reached or canceled.
      for (ResultsIterations = 1; ResultsIterations < MaximumGenerations; ResultsIterations++) {
        try {
          Iterate();
          cancellationToken.ThrowIfCancellationRequested();
        }
        finally {
          AnalyzeSolutions();
        }
      }
    }
    private void Iterate() {
      var offspring = MutateOffspring();
      PenalizeEvaluate(offspring);
      var parents = solutions.Concat(offspring).ToArray();
      SelectParents(parents, solutions.Length);
      UpdatePopulation(parents);
    }
    protected override void OnExecutionTimeChanged() {
      base.OnExecutionTimeChanged();
      if (CancellationTokenSource == null) return;
      if (MaximumRuntime == -1) return;
      if (ExecutionTime.TotalSeconds > MaximumRuntime) CancellationTokenSource.Cancel();
    }
    #endregion

    #region Initialization
    private MOCMAESIndividual InitializeIndividual(RealVector x) {
      var zeros = new RealVector(x.Length);
      var identity = new double[x.Length, x.Length];
      for (var i = 0; i < x.Length; i++) identity[i, i] = 1;
      return new MOCMAESIndividual(x, internals.TargetSuccessProbability, InitialSigma, zeros, identity, internals);
    }
    private void InitSolutions() {
      solutions = new MOCMAESIndividual[PopulationSize];
      for (var i = 0; i < PopulationSize; i++) {
        var x = new RealVector(Problem.ProblemSize); // Uniform distibution in all dimensions assumed.
                                                     // There is the UniformSolutionCreater associated with the Encoding, but it was considered not usable here
        var bounds = Problem.Bounds;
        for (var j = 0; j < Problem.ProblemSize; j++) {
          var dim = j % bounds.Rows;
          x[j] = random.NextDouble() * (bounds[dim, 1] - bounds[dim, 0]) + bounds[dim, 0];
        }
        solutions[i] = InitializeIndividual(x);
      }
      PenalizeEvaluate(solutions);
    }
    private void InitStrategy() {
      const int lambda = 1;
      double n = Problem.ProblemSize;

      var targetSuccessProbability = 1.0 / (5.0 + Math.Sqrt(lambda) / 2.0);
      var stepSizeDampeningFactor = 1.0 + n / (2.0 * lambda);
      var stepSizeLearningRate = targetSuccessProbability * lambda / (2.0 + targetSuccessProbability * lambda);
      var evolutionPathLearningRate = 2.0 / (n + 2.0);
      var covarianceMatrixLearningRate = 2.0 / (n * n + 6.0);
      var covarianceMatrixUnlearningRate = 0.4 / (Math.Pow(n, 1.6) + 1);
      var successThreshold = 0.44;
      internals = new MOCMAESParameters(stepSizeLearningRate, stepSizeDampeningFactor, targetSuccessProbability, evolutionPathLearningRate, covarianceMatrixLearningRate, covarianceMatrixUnlearningRate, successThreshold);

    }
    private void InitResults() {
      AddResult(IterationsResultName, "The number of gererations evaluated", new IntValue(0));
      AddResult(EvaluationsResultName, "The number of function evaltions performed", new IntValue(0));
      AddResult(HypervolumeResultName, "The hypervolume of the current front considering the Referencepoint defined in the Problem", new DoubleValue(0.0));
      AddResult(BestHypervolumeResultName, "The best hypervolume of the current run considering the Referencepoint defined in the Problem", new DoubleValue(0.0));
      AddResult(BestKnownHypervolumeResultName, "The best knwon hypervolume considering the Referencepoint defined in the Problem", new DoubleValue(double.NaN));
      AddResult(DifferenceToBestKnownHypervolumeResultName, "The difference between the current and the best known hypervolume", new DoubleValue(double.NaN));
      AddResult(GenerationalDistanceResultName, "The generational distance to an optimal pareto front defined in the Problem", new DoubleValue(double.NaN));
      AddResult(InvertedGenerationalDistanceResultName, "The inverted generational distance to an optimal pareto front defined in the Problem", new DoubleValue(double.NaN));
      AddResult(CrowdingResultName, "The average crowding value for the current front (excluding infinities)", new DoubleValue(0.0));
      AddResult(SpacingResultName, "The spacing for the current front (excluding infinities)", new DoubleValue(0.0));

      var table = new DataTable("QualityIndicators");
      table.Rows.Add(new DataRow(BestHypervolumeResultName));
      table.Rows.Add(new DataRow(HypervolumeResultName));
      table.Rows.Add(new DataRow(CrowdingResultName));
      table.Rows.Add(new DataRow(GenerationalDistanceResultName));
      table.Rows.Add(new DataRow(InvertedGenerationalDistanceResultName));
      table.Rows.Add(new DataRow(DifferenceToBestKnownHypervolumeResultName));
      table.Rows.Add(new DataRow(SpacingResultName));
      AddResult(TimetableResultName, "Different quality meassures in a timeseries", table);
      AddResult(SolutionsResultName, "The current front", new DoubleMatrix());
      AddResult(ScatterPlotResultName, "A scatterplot displaying the evaluated solutions and (if available) the analytically optimal front", new ScatterPlotContent(null, null, null, 2));

      if (Problem.BestKnownFront != null) {
        ResultsBestKnownHypervolume = Hypervolume.Calculate(Utilities.ToArray(Problem.BestKnownFront), Problem.TestFunction.ReferencePoint(Problem.Objectives), Problem.Maximization);
        ResultsDifferenceBestKnownHypervolume = ResultsBestKnownHypervolume;
      }
      ResultsScatterPlot = new ScatterPlotContent(null, null, Utilities.ToArray(Problem.BestKnownFront), Problem.Objectives);
    }
    #endregion

    private MOCMAESIndividual[] MutateOffspring() {
      var offspring = new MOCMAESIndividual[PopulationSize];
      for (var i = 0; i < PopulationSize; i++) {
        offspring[i] = new MOCMAESIndividual(solutions[i]);
        offspring[i].Mutate(gauss);
      }
      return offspring;
    }

    #region Evaluation
    private void PenalizeEvaluate(IEnumerable<MOCMAESIndividual> offspring) {
      foreach (var child in offspring) {
        if (IsFeasable(child.Mean)) {
          child.Fitness = Evaluate(child.Mean);
          child.PenalizedFitness = child.Fitness;
        } else {
          var t = ClosestFeasible(child.Mean);
          child.Fitness = Evaluate(t);
          child.PenalizedFitness = Penalize(child.Mean, t, child.Fitness);
        }
      }
    }
    private double[] Evaluate(RealVector x) {
      var res = Problem.Evaluate(x);
      ResultsEvaluations++;
      return res;
    }
    private static double[] Penalize(RealVector x, RealVector t, IEnumerable<double> fitness) {
      var penalty = x.Select((t1, i) => t1 - t[i]).Sum(d => d * d);
      return fitness.Select(v => v + penalty).ToArray();
    }
    private RealVector ClosestFeasible(RealVector x) {
      var bounds = Problem.Bounds;
      var r = new RealVector(x.Length);
      for (var i = 0; i < x.Length; i++) {
        var dim = i % bounds.Rows;
        r[i] = Math.Min(Math.Max(bounds[dim, 0], x[i]), bounds[dim, 1]);
      }
      return r;
    }
    private bool IsFeasable(RealVector offspring) {
      var bounds = Problem.Bounds;
      for (var i = 0; i < offspring.Length; i++) {
        var dim = i % bounds.Rows;
        if (bounds[dim, 0] > offspring[i] || offspring[i] > bounds[dim, 1]) return false;
      }
      return true;
    }
    #endregion

    private void SelectParents(IReadOnlyList<MOCMAESIndividual> parents, int length) {
      //perform a nondominated sort to assign the rank to every element
      var fronts = NonDominatedSort(parents);

      //deselect the highest rank fronts until we would end up with less or equal mu elements
      var rank = fronts.Count - 1;
      var popSize = parents.Count;
      while (popSize - fronts[rank].Count >= length) {
        var front = fronts[rank];
        foreach (var i in front) i.Selected = false;
        popSize -= front.Count;
        rank--;
      }

      //now use the indicator to deselect the approximatingly worst elements of the last selected front
      var front1 = fronts[rank];
      for (; popSize > length; popSize--) {
        var lc = Indicator.LeastContributer(front1.ToArray(), x => x.PenalizedFitness, Problem);
        front1[lc].Selected = false;
        front1.Swap(lc, front1.Count - 1);
        front1.RemoveAt(front1.Count - 1);
      }
    }

    private void UpdatePopulation(IReadOnlyList<MOCMAESIndividual> parents) {
      var offspringSucess = new int[solutions.Length];
      var offspringLength = parents.Count - solutions.Length;
      for (var i = 0; i < offspringLength; i++) {
        if (!parents[i + solutions.Length].Selected) continue;
        parents[i + solutions.Length].UpdateAsOffspring();
        offspringSucess[i] += MOCMAESIndividual.Success;
      }
      for (var i = 0; i < solutions.Length; i++) if (parents[i].Selected) parents[i].UpdateAsParent(offspringSucess[i]);
      solutions = new MOCMAESIndividual[solutions.Length];
      var j = 0;
      foreach (var ind in parents) if (ind.Selected) solutions[j++] = ind;
    }

    #region Analysis
    private void AnalyzeSolutions() {
      ResultsScatterPlot = new ScatterPlotContent(solutions.Select(x => x.Fitness).ToArray(),
          solutions.Select(x => x.Mean.ToArray()).ToArray(),
          ResultsScatterPlot.ParetoFront,
          ResultsScatterPlot.Objectives);
      ResultsSolutions = ToMatrix(solutions.Select(x => x.Mean.ToArray()));
      AnalyzeQualityIndicators();
    }
    private void AnalyzeQualityIndicators() {
      var front = NonDominatedSelect.GetDominatingVectors(solutions.Select(x => x.Fitness), Problem.TestFunction.ReferencePoint(Problem.Objectives), Problem.Maximization, true).ToArray();
      var bounds = Problem.Bounds.CloneAsMatrix();
      ResultsCrowding = Crowding.Calculate(front, bounds);
      ResultsSpacing = Spacing.Calculate(front);
      ResultsGenerationalDistance = Problem.BestKnownFront != null ? GenerationalDistance.Calculate(front, Utilities.ToArray(Problem.BestKnownFront), 1) : double.NaN;

      ResultsInvertedGenerationalDistance = Problem.BestKnownFront != null ? InvertedGenerationalDistance.Calculate(front, Utilities.ToArray(Problem.BestKnownFront), 1) : double.NaN;

      ResultsHypervolume = Hypervolume.Calculate(front, Problem.TestFunction.ReferencePoint(Problem.Objectives), Problem.Maximization);
      ResultsBestHypervolume = Math.Max(ResultsHypervolume, ResultsBestHypervolume);
      ResultsDifferenceBestKnownHypervolume = ResultsBestKnownHypervolume - ResultsBestHypervolume;

      //Datalines
      ResultsBestHypervolumeDataLine.Values.Add(ResultsBestHypervolume);
      ResultsHypervolumeDataLine.Values.Add(ResultsHypervolume);
      ResultsCrowdingDataLine.Values.Add(ResultsCrowding);
      ResultsGenerationalDistanceDataLine.Values.Add(ResultsGenerationalDistance);
      ResultsInvertedGenerationalDistanceDataLine.Values.Add(ResultsInvertedGenerationalDistance);
      ResultsSpacingDataLine.Values.Add(ResultsSpacing);
      ResultsHypervolumeDifferenceDataLine.Values.Add(ResultsDifferenceBestKnownHypervolume);

    }
    #endregion



    #region Helpers
    public DoubleMatrix ToMatrix(IEnumerable<double[]> data) {
      var d2 = data.ToArray();
      var mat = new DoubleMatrix(d2.Length, d2[0].Length);
      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 void AddResult<T>(string name1, string desc, T defaultValue) where T : class, IItem {
      Results.Add(new Result(name1, desc, defaultValue));
    }
    //blatantly stolen form HeuristicLab.Optimization.Operators.FastNonDominatedSort
    //however: Operators.FastNonDominatedSort does not return ranked fronts => rerank after sorting would not be sensible
    #region FastNonDominatedSort   
    private enum DominationResult { Dominates, IsDominated, IsNonDominated };
    private List<List<MOCMAESIndividual>> NonDominatedSort(IReadOnlyList<MOCMAESIndividual> individuals) {
      const bool dominateOnEqualQualities = false;
      var maximization = Problem.Maximization;
      if (individuals == null) throw new InvalidOperationException(Name + ": No qualities found.");
      var populationSize = individuals.Count;

      var fronts = new List<List<MOCMAESIndividual>>();
      var dominatedScopes = new Dictionary<MOCMAESIndividual, List<int>>();
      var dominationCounter = new int[populationSize];

      for (var pI = 0; pI < populationSize - 1; pI++) {
        var p = individuals[pI];
        List<int> dominatedScopesByp;
        if (!dominatedScopes.TryGetValue(p, out dominatedScopesByp))
          dominatedScopes[p] = dominatedScopesByp = new List<int>();
        for (var qI = pI + 1; qI < populationSize; qI++) {
          var test = Dominates(individuals[pI], individuals[qI], maximization, dominateOnEqualQualities);
          if (test == DominationResult.Dominates) {
            dominatedScopesByp.Add(qI);
            dominationCounter[qI] += 1;
          } else if (test == DominationResult.IsDominated) {
            dominationCounter[pI] += 1;
            if (!dominatedScopes.ContainsKey(individuals[qI]))
              dominatedScopes.Add(individuals[qI], new List<int>());
            dominatedScopes[individuals[qI]].Add(pI);
          }
          if (pI == populationSize - 2
            && qI == populationSize - 1
            && dominationCounter[qI] == 0) {
            AddToFront(individuals[qI], fronts, 0);
          }
        }
        if (dominationCounter[pI] == 0) {
          AddToFront(p, fronts, 0);
        }
      }
      var i = 0;
      while (i < fronts.Count && fronts[i].Count > 0) {
        var nextFront = new List<MOCMAESIndividual>();
        foreach (var p in fronts[i]) {
          List<int> dominatedScopesByp;
          if (!dominatedScopes.TryGetValue(p, out dominatedScopesByp)) continue;
          foreach (var dominatedScope in dominatedScopesByp) {
            dominationCounter[dominatedScope] -= 1;
            if (dominationCounter[dominatedScope] != 0) continue;
            nextFront.Add(individuals[dominatedScope]);
          }
        }
        i += 1;
        fronts.Add(nextFront);
      }

      for (i = 0; i < fronts.Count; i++) {
        foreach (var p in fronts[i]) {
          p.Rank = i;
        }
      }
      return fronts;
    }
    private static void AddToFront(MOCMAESIndividual p, IList<List<MOCMAESIndividual>> fronts, int i) {
      if (i == fronts.Count) fronts.Add(new List<MOCMAESIndividual>());
      fronts[i].Add(p);
    }
    private static DominationResult Dominates(MOCMAESIndividual left, MOCMAESIndividual right, bool[] maximizations, bool dominateOnEqualQualities) {
      return Dominates(left.PenalizedFitness, right.PenalizedFitness, maximizations, dominateOnEqualQualities);
    }
    private static DominationResult Dominates(IReadOnlyList<double> left, IReadOnlyList<double> right, IReadOnlyList<bool> maximizations, bool dominateOnEqualQualities) {
      //mkommend Caution: do not use LINQ.SequenceEqual for comparing the two quality arrays (left and right) due to performance reasons
      if (dominateOnEqualQualities) {
        var equal = true;
        for (var i = 0; i < left.Count; i++) {
          if (left[i] != right[i]) {
            equal = false;
            break;
          }
        }
        if (equal) return DominationResult.Dominates;
      }

      bool leftIsBetter = false, rightIsBetter = false;
      for (var i = 0; i < left.Count; i++) {
        if (IsDominated(left[i], right[i], maximizations[i])) rightIsBetter = true;
        else if (IsDominated(right[i], left[i], maximizations[i])) leftIsBetter = true;
        if (leftIsBetter && rightIsBetter) break;
      }

      if (leftIsBetter && !rightIsBetter) return DominationResult.Dominates;
      if (!leftIsBetter && rightIsBetter) return DominationResult.IsDominated;
      return DominationResult.IsNonDominated;
    }
    private static bool IsDominated(double left, double right, bool maximization) {
      return maximization && left < right
        || !maximization && left > right;
    }
    #endregion
    #endregion
  }
}