source: branches/HeuristicLab.Algorithms.CMAEvolutionStrategy/sources/3.3/MOCMAES/MOCMASEvolutionStrategy.cs @ 13793

#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.Threading;
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.MultiObjectiveTestFunctions;
using HeuristicLab.Random;

namespace HeuristicLab.Algorithms.CMAEvolutionStrategy {
  [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; }
    }

    private readonly IRandom random = new MersenneTwister();


    #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 InitialIterationsName = "InitialIterations";
    private const string InitialSigmaName = "InitialSigma";
    private const string MuName = "Mu";
    private const string CMAInitializerName = "CMAInitializer";
    private const string CMAMutatorName = "CMAMutator";
    private const string CMARecombinatorName = "CMARecombinator";
    private const string CMAUpdaterName = "CMAUpdater";
    private const string AnalyzerName = "Analyzer";
    private const string MaximumGenerationsName = "MaximumGenerations";
    private const string MaximumEvaluatedSolutionsName = "MaximumEvaluatedSolutions";
    private const string TargetQualityName = "TargetQuality";
    private const string MinimumQualityChangeName = "MinimumQualityChange";
    private const string MinimumQualityHistoryChangeName = "MinimumQualityHistoryChange";
    private const string MinimumStandardDeviationName = "MinimumStandardDeviation";
    private const string MaximumStandardDeviationChangeName = "MaximumStandardDeviationChange";
    #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]; }
    }
    public IValueParameter<MultiAnalyzer> AnalyzerParameter {
      get { return (IValueParameter<MultiAnalyzer>)Parameters[AnalyzerName]; }
    }
    private IFixedValueParameter<IntValue> InitialIterationsParameter {
      get { return (IFixedValueParameter<IntValue>)Parameters[InitialIterationsName]; }
    }
    public IValueParameter<DoubleValue> InitialSigmaParameter {
      get { return (IValueParameter<DoubleValue>)Parameters[InitialSigmaName]; }
    }
    private OptionalValueParameter<IntValue> MuParameter {
      get { return (OptionalValueParameter<IntValue>)Parameters[MuName]; }
    }
    public IConstrainedValueParameter<ICMAInitializer> CMAInitializerParameter {
      get { return (IConstrainedValueParameter<ICMAInitializer>)Parameters[CMAInitializerName]; }
    }
    public IConstrainedValueParameter<ICMAManipulator> CMAMutatorParameter {
      get { return (IConstrainedValueParameter<ICMAManipulator>)Parameters[CMAMutatorName]; }
    }
    public IConstrainedValueParameter<ICMARecombinator> CMARecombinatorParameter {
      get { return (IConstrainedValueParameter<ICMARecombinator>)Parameters[CMARecombinatorName]; }
    }
    public IConstrainedValueParameter<ICMAUpdater> CMAUpdaterParameter {
      get { return (IConstrainedValueParameter<ICMAUpdater>)Parameters[CMAUpdaterName]; }
    }
    private IFixedValueParameter<IntValue> MaximumGenerationsParameter {
      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumGenerationsName]; }
    }
    private IFixedValueParameter<IntValue> MaximumEvaluatedSolutionsParameter {
      get { return (IFixedValueParameter<IntValue>)Parameters[MaximumEvaluatedSolutionsName]; }
    }
    private IFixedValueParameter<DoubleValue> TargetQualityParameter {
      get { return (IFixedValueParameter<DoubleValue>)Parameters[TargetQualityName]; }
    }
    private IFixedValueParameter<DoubleValue> MinimumQualityChangeParameter {
      get { return (IFixedValueParameter<DoubleValue>)Parameters[MinimumQualityChangeName]; }
    }
    private IFixedValueParameter<DoubleValue> MinimumQualityHistoryChangeParameter {
      get { return (IFixedValueParameter<DoubleValue>)Parameters[MinimumQualityHistoryChangeName]; }
    }
    private IFixedValueParameter<DoubleValue> MinimumStandardDeviationParameter {
      get { return (IFixedValueParameter<DoubleValue>)Parameters[MinimumStandardDeviationName]; }
    }
    private IFixedValueParameter<DoubleValue> MaximumStandardDeviationChangeParameter {
      get { return (IFixedValueParameter<DoubleValue>)Parameters[MaximumStandardDeviationChangeName]; }
    }
    #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 InitialIterations {
      get { return InitialIterationsParameter.Value.Value; }
      set { InitialIterationsParameter.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 TargetQuality {
      get { return TargetQualityParameter.Value.Value; }
      set { TargetQualityParameter.Value.Value = value; }
    }
    public double MinimumQualityChange {
      get { return MinimumQualityChangeParameter.Value.Value; }
      set { MinimumQualityChangeParameter.Value.Value = value; }
    }
    public double MinimumQualityHistoryChange {
      get { return MinimumQualityHistoryChangeParameter.Value.Value; }
      set { MinimumQualityHistoryChangeParameter.Value.Value = value; }
    }
    public double MinimumStandardDeviation {
      get { return MinimumStandardDeviationParameter.Value.Value; }
      set { MinimumStandardDeviationParameter.Value.Value = value; }
    }
    public double MaximumStandardDeviationChange {
      get { return MaximumStandardDeviationChangeParameter.Value.Value; }
      set { MaximumStandardDeviationChangeParameter.Value.Value = value; }
    }
    public DoubleValue InitialSigma {
      get { return InitialSigmaParameter.Value; }
      set { InitialSigmaParameter.Value = value; }
    }
    public IntValue Mu {
      get { return MuParameter.Value; }
      set { MuParameter.Value = value; }
    }
    public ICMAInitializer CMAInitializer {
      get { return CMAInitializerParameter.Value; }
      set { CMAInitializerParameter.Value = value; }
    }
    public ICMAManipulator CMAMutator {
      get { return CMAMutatorParameter.Value; }
      set { CMAMutatorParameter.Value = value; }
    }
    public ICMARecombinator CMARecombinator {
      get { return CMARecombinatorParameter.Value; }
      set { CMARecombinatorParameter.Value = value; }
    }
    public MultiAnalyzer Analyzer {
      get { return AnalyzerParameter.Value; }
      set { AnalyzerParameter.Value = value; }
    }
    public ICMAUpdater CMAUpdater {
      get { return CMAUpdaterParameter.Value; }
      set { CMAUpdaterParameter.Value = value; }
    }

    #endregion

    #region ResultsProperties
    private double ResultsBestQuality {
      get { return ((DoubleValue)Results["Best Quality"].Value).Value; }
      set { ((DoubleValue)Results["Best Quality"].Value).Value = value; }
    }

    private RealVector ResultsBestSolution {
      get { return (RealVector)Results["Best Solution"].Value; }
      set { Results["Best Solution"].Value = value; }
    }

    private int ResultsBestFoundOnEvaluation {
      get { return ((IntValue)Results["Evaluation Best Solution Was Found"].Value).Value; }
      set { ((IntValue)Results["Evaluation Best Solution Was Found"].Value).Value = value; }
    }

    private int ResultsEvaluations {
      get { return ((IntValue)Results["Evaluations"].Value).Value; }
      set { ((IntValue)Results["Evaluations"].Value).Value = value; }
    }
    private int ResultsIterations {
      get { return ((IntValue)Results["Iterations"].Value).Value; }
      set { ((IntValue)Results["Iterations"].Value).Value = value; }
    }

    private DataTable ResultsQualities {
      get { return ((DataTable)Results["Qualities"].Value); }
    }
    private DataRow ResultsQualitiesBest {
      get { return ResultsQualities.Rows["Best Quality"]; }
    }

    private DataRow ResultsQualitiesIteration {
      get { return ResultsQualities.Rows["Iteration Quality"]; }
    }

    #endregion

    #region constructors and hlBoilerPlate-code
    [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);
    }

    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<IntValue>(InitialIterationsName, "The number of iterations that should be performed with only axis parallel mutation.", new IntValue(0)));
      Parameters.Add(new FixedValueParameter<DoubleValue>(InitialSigmaName, "The initial sigma can be a single value or a value for each dimension. All values need to be > 0.", new DoubleValue(0.5)));
      Parameters.Add(new OptionalValueParameter<IntValue>(MuName, "Optional, the mu best offspring that should be considered for update of the new mean and strategy parameters. If not given it will be automatically calculated."));
      Parameters.Add(new ConstrainedValueParameter<ICMARecombinator>(CMARecombinatorName, "The operator used to calculate the new mean."));
      Parameters.Add(new ConstrainedValueParameter<ICMAManipulator>(CMAMutatorName, "The operator used to manipulate a point."));
      Parameters.Add(new ConstrainedValueParameter<ICMAInitializer>(CMAInitializerName, "The operator that initializes the covariance matrix and strategy parameters."));
      Parameters.Add(new ConstrainedValueParameter<ICMAUpdater>(CMAUpdaterName, "The operator that updates the covariance matrix and strategy parameters."));
      Parameters.Add(new ValueParameter<MultiAnalyzer>(AnalyzerName, "The operator used to analyze each generation.", new MultiAnalyzer()));
      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)));
      Parameters.Add(new FixedValueParameter<DoubleValue>(TargetQualityName, "(stopFitness) Surpassing this quality value terminates the algorithm.", new DoubleValue(double.NaN)));
      Parameters.Add(new FixedValueParameter<DoubleValue>(MinimumQualityChangeName, "(stopTolFun) If the range of fitness values is less than a certain value the algorithm terminates (set to 0 or positive value to enable).", new DoubleValue(double.NaN)));
      Parameters.Add(new FixedValueParameter<DoubleValue>(MinimumQualityHistoryChangeName, "(stopTolFunHist) If the range of fitness values is less than a certain value for a certain time the algorithm terminates (set to 0 or positive to enable).", new DoubleValue(double.NaN)));
      Parameters.Add(new FixedValueParameter<DoubleValue>(MinimumStandardDeviationName, "(stopTolXFactor) If the standard deviation falls below a certain value the algorithm terminates (set to 0 or positive to enable).", new DoubleValue(double.NaN)));
      Parameters.Add(new FixedValueParameter<DoubleValue>(MaximumStandardDeviationChangeName, "(stopTolUpXFactor) If the standard deviation changes by a value larger than this parameter the algorithm stops (set to a value > 0 to enable).", new DoubleValue(double.NaN)));

    }
    #endregion

    #region updates 
    private void updatePopulation(CMAHansenIndividual[] parents) {
      int[] offspringSucess = new int[solutions.Length];
      int offspringLength = parents.Length - solutions.Length;

      for (int i = 0; i < offspringLength; i++) {
        if (parents[i + solutions.Length].selected) {
          updateAsOffspring(parents[i + solutions.Length]);

          //TODO this may change if more offspring per parent is allowed 
          offspringSucess[i] += CMAHansenIndividual.success;
        }
      }
      for (int i = 0; i < solutions.Length; i++) {
        if (parents[i].selected) {
          updateAsParent(parents[i], offspringSucess[i]);
        }
      }

      solutions = new CMAHansenIndividual[solutions.Length];
      int j = 0;
      foreach (CMAHansenIndividual ind in parents) {
        if (ind.selected) {
          solutions[j++] = ind;
        }
      }

    }

    private void updateAsParent(CMAHansenIndividual c, int v) {
      c.successProbability = (1 - stepSizeLearningRate) * c.successProbability + stepSizeLearningRate * (v == CMAHansenIndividual.success ? 1 : 0);
      c.sigma *= Math.Exp(1 / stepSizeDampeningFactor * (c.successProbability - targetSuccessProbability) / (1 - targetSuccessProbability));
      if (v != CMAHansenIndividual.failure) return;
      if (c.successProbability < successThreshold) {
        double stepNormSqr = c.getSetpNormSqr();
        double rate = covarianceMatrixUnlearningRate;
        if (stepNormSqr > 1 && 1 < covarianceMatrixUnlearningRate * (2 * stepNormSqr - 1)) {
          rate = 1 / (2 * stepNormSqr - 1);
        }
        rankOneUpdate(c, 1 + rate, -rate, c.lastStep);

      } else {
        roundUpdate(c);
      }

    }

    private void updateAsOffspring(CMAHansenIndividual c) {
      c.successProbability = (1 - stepSizeLearningRate) * c.successProbability + stepSizeLearningRate;
      c.sigma *= Math.Exp(1 / stepSizeDampeningFactor * (c.successProbability - targetSuccessProbability) / (1 - targetSuccessProbability));
      double evolutionpathUpdateWeight = evolutionPathLearningRate * (2.0 - evolutionPathLearningRate);
      if (c.successProbability < successThreshold) {
        c.updateEvolutionPath(1 - evolutionPathLearningRate, evolutionpathUpdateWeight);
        rankOneUpdate(c, 1 - covarianceMatrixLearningRate, covarianceMatrixLearningRate, c.evolutionPath);
      } else {
        roundUpdate(c);
      }
    }

    private void rankOneUpdate(CMAHansenIndividual c, double v1, double v2, RealVector lastStep) {
      c.choleskyUpdate(lastStep, v1, v2);
    }

    private void roundUpdate(CMAHansenIndividual c) {
      double evolutionPathUpdateWeight = evolutionPathLearningRate * (2.0 - evolutionPathLearningRate);
      c.updateEvolutionPath(1 - evolutionPathLearningRate, 0);
      rankOneUpdate(c, 1 - covarianceMatrixLearningRate + evolutionPathUpdateWeight,
        covarianceMatrixLearningRate, c.evolutionPath);
    }
    #endregion

    #region selection

    private T[] merge<T>(T[] parents, T[] offspring) {
      T[] merged = new T[parents.Length + offspring.Length];
      for (int i = 0; i < parents.Length; i++) {
        merged[i] = parents[i];
      }
      for (int i = 0; i < offspring.Length; i++) {
        merged[i + parents.Length] = offspring[i];
      }
      return merged;
    }

    private void selection(CMAHansenIndividual[] 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
      int rank = fronts.Count - 1;
      int popSize = parents.Length;
      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 worst approximating elements of the last selected front
      var front_ = fronts[rank];
      for (; popSize > length; popSize--) {
        int lc = indicator.leastContributer<CMAHansenIndividual>(front_.ToArray(), x => x.penalizedFitness);   //TODO: This is a battel in its own right to be fought another day
        front_[lc].selected = false;
        front_.Swap(lc, front_.Count - 1);
        fronts.RemoveAt(front_.Count - 1);
      }
    }
    #endregion

    #region penalize Box-Constraints
    private void penalizeEvaluate(IEnumerable<CMAHansenIndividual> offspring) {
      foreach (CMAHansenIndividual child in offspring) {
        penalizeEvaluate(child);
      }
    }

    private void penalizeEvaluate(CMAHansenIndividual offspring) {
      if (isFeasable(offspring.x)) {
        offspring.fitness = Problem.Evaluate(offspring.x, random);
        offspring.penalizedFitness = offspring.fitness;
      } else {
        RealVector t = closestFeasible(offspring.x);
        offspring.fitness = Problem.Evaluate(t, random);
        offspring.penalizedFitness = penalize(offspring.x, t, offspring.fitness);
      }
    }

    private double[] penalize(RealVector x, RealVector t, double[] fitness) {
      double[] d = new double[fitness.Length];
      double penality = penalize(x, t);
      for (int i = 0; i < fitness.Length; i++) {
        d[i] = fitness[i] + penality;
      }
      return d;
    }

    private double penalize(RealVector x, RealVector t) {
      double sum = 0;
      for (int i = 0; i < x.Length; i++) {
        double d = x[i] - t[i];
        sum += d * d;
      }
      return sum;
    }

    private RealVector closestFeasible(RealVector x) {
      DoubleMatrix bounds = Problem.Bounds;
      RealVector r = new RealVector(x.Length);
      for (int i = 0; i < x.Length; i++) {
        r[i] = Math.Min(Math.Max(bounds[i, 0], x[i]), bounds[i, 1]);
      }
      return r;
    }

    private bool isFeasable(RealVector offspring) {
      DoubleMatrix bounds = Problem.Bounds;
      for (int i = 0; i < offspring.Length; i++) {
        if (bounds[i, 0] > offspring[i] || offspring[i] > bounds[i, 1]) return false;
      }
      return true;
    }

    #endregion

    #region mutation
    private CMAHansenIndividual[] generateOffspring() {
      CMAHansenIndividual[] offspring = new CMAHansenIndividual[PopulationSize];  //TODO this changes if 1,1-ES is replaced with 1,n-ES
      for (int i = 0; i < offspring.Length; i++) {
        offspring[i] = new CMAHansenIndividual(solutions[i]);
        offspring[i].mutate(gauss);
      }
      return offspring;
    }
    #endregion

    #region initialization
    private CMAHansenIndividual initializeIndividual(RealVector x) {
      var zeros = new RealVector(x.Length);
      var identity = new double[x.Length, x.Length];
      for (int i = 0; i < x.Length; i++) {
        identity[i, i] = 1;
      }
      return new CMAHansenIndividual(x, targetSuccessProbability, InitialSigma.Value, zeros, identity);
    }

    private void initSolutions() {
      for (int i = 0; i < PopulationSize; i++) {
        RealVector x = null; //TODO get those with magic 
        solutions[i] = initializeIndividual(x);
      }
      penalizeEvaluate(solutions);
    }

    private void initStrategy() {
      int lambda = 1;
      double n = Problem.ProblemSize;
      indicator = new CrowdingIndicator();
      gauss = new NormalDistributedRandom(random, 0, 1);
      targetSuccessProbability = 1.0 / (5.0 + Math.Sqrt(lambda) / 2.0);
      stepSizeDampeningFactor = 1.0 + n / (2.0 * lambda);
      stepSizeLearningRate = targetSuccessProbability * lambda / (2.0 + targetSuccessProbability * lambda);
      evolutionPathLearningRate = 2.0 / (n + 2.0);
      covarianceMatrixLearningRate = 2.0 / (n * n + 6.0);
      covarianceMatrixUnlearningRate = 0.4 / (Math.Pow(n, 1.6) + 1);
      successThreshold = 0.44;

    }

    #endregion

    #region analyze
    private void analyzeSolutions() {
      //TODO give HyperVolume stuff
      //Problem.Analyze();   

    }

    #endregion


    //MU = populationSize
    #region mainloop
    protected override void Run(CancellationToken cancellationToken) {
      // Set up the algorithm
      if (SetSeedRandomly) Seed = new System.Random().Next();
      random.Reset(Seed);


      // Set up the results display
      Results.Add(new Result("Iterations", new IntValue(0)));
      Results.Add(new Result("Evaluations", new IntValue(0)));
      var table = new DataTable("Hypervolumes");
      table.Rows.Add(new DataRow("Best Hypervolumes"));
      var iterationRows = new DataRow("Iteration Hypervolumes");
      iterationRows.VisualProperties.LineStyle = DataRowVisualProperties.DataRowLineStyle.Dot;
      table.Rows.Add(iterationRows);
      Results.Add(new Result("Hypervolumes", table));


      initStrategy();
      initSolutions();



      // Loop until iteration limit reached or canceled.
      for (ResultsIterations = 0; ResultsIterations < MaximumGenerations; ResultsIterations++) {

        try {
          iterate();
          cancellationToken.ThrowIfCancellationRequested();
        }
        finally {
        }
      }
    }

    protected override void OnExecutionTimeChanged() {
      base.OnExecutionTimeChanged();
      if (CancellationTokenSource == null) return;
      if (MaximumRuntime == -1) return;
      if (ExecutionTime.TotalSeconds > MaximumRuntime) CancellationTokenSource.Cancel();
    }

    private void iterate() {
      CMAHansenIndividual[] offspring = generateOffspring();
      penalizeEvaluate(offspring);
      selection(merge(solutions, offspring), solutions.Length);
      updatePopulation(offspring);

    }
    #endregion

    #region bernhard properties
    private NormalDistributedRandom gauss;
    private CMAHansenIndividual[] solutions;
    private const double penalizeFactor = 1e-6;
    private IIndicator indicator;

    private double stepSizeLearningRate; //=cp learning rate in [0,1]
    private double stepSizeDampeningFactor; //d
    private double targetSuccessProbability;// p^target_succ
    private double evolutionPathLearningRate;//cc
    private double covarianceMatrixLearningRate;//ccov
    private double covarianceMatrixUnlearningRate;   //from shark
    private double successThreshold; //ptresh
    #endregion

    private class CMAHansenIndividual {
      public static readonly int success = 1;
      public static readonly int noSuccess = 2;
      public static readonly int failure = 3;


      //Chromosome
      public RealVector x;
      public double successProbability;
      public double sigma;//stepsize
      public RealVector evolutionPath; // pc
      public RealVector lastStep;
      public RealVector lastZ;
      private double[,] lowerCholesky;


      //Phenotype
      public double[] fitness;
      public double[] penalizedFitness;
      public bool selected = true;

      internal double rank;



      /// <summary>
      /// 
      /// </summary>
      /// <param name="x">has to be 0-vector with correct lenght</param>
      /// <param name="p_succ">has to be ptargetsucc</param>
      /// <param name="sigma">initialSigma</param>
      /// <param name="pc">has to be 0-vector with correct lenght</param>
      /// <param name="C">has to be a symmetric positive definit Covariance matrix</param>
      public CMAHansenIndividual(RealVector x, double p_succ, double sigma, RealVector pc, double[,] C) {
        this.x = x;
        this.successProbability = p_succ;
        this.sigma = sigma;
        this.evolutionPath = pc;
        choleskyDecomposition(C);
      }

      private void choleskyDecomposition(double[,] C) {
        if (C.GetLength(0) != C.GetLength(1)) throw new ArgumentException("Covariancematrix is not quadratic");
        int n = C.GetLength(0);
        lowerCholesky = new double[n, n];
        double[,] A = lowerCholesky;
        for (int i = 0; i < n; i++) {
          for (int j = 0; j <= i; j++) {
            A[i, j] = C[i, j];   //simulate inplace transform
            double sum = A[i, j];
            for (int k = 0; k < j; k++) {
              sum = sum - A[i, k] * A[j, k];
              if (i > j) { A[i, j] = sum / A[j, j]; } else if (sum > 0) {
                A[i, i] = Math.Sqrt(sum);
              } else {
                throw new ArgumentException("Covariancematrix is not symetric positiv definit");
              }
            }
          }

        }
      }

      public CMAHansenIndividual(CMAHansenIndividual other) {


        this.successProbability = other.successProbability;
        this.sigma = other.sigma;

        // no no no make DEEP copies
        this.x = (RealVector)other.x.Clone();    //This may be ommited
        this.evolutionPath = (RealVector)other.evolutionPath.Clone();
        this.lastStep = (RealVector)lastStep.Clone();  //This may be ommited
        this.lastZ = (RealVector)lastZ.Clone();  //This may be ommited


        this.lowerCholesky = (double[,])other.lowerCholesky.Clone();
      }

      public void updateEvolutionPath(double learningRate, double updateWeight) {
        updateWeight = Math.Sqrt(updateWeight);
        for (int i = 0; i < evolutionPath.Length; i++) {
          evolutionPath[i] *= learningRate;
          evolutionPath[i] += updateWeight * lastStep[i];
        }
      }

      public double getSetpNormSqr() {
        double sum = 0;
        foreach (double d in lastZ) {
          sum += d * d;
        }
        return sum;
      }

      public void choleskyUpdate(RealVector v, double alpha, double beta) {
        int n = v.Length;
        double[] temp = new double[n];
        for (int i = 0; i < n; i++) temp[i] = v[i];
        double betaPrime = 1;
        double a = Math.Sqrt(alpha);
        for (int j = 0; j < n; j++) {
          double Ljj = a * lowerCholesky[j, j];
          double dj = Ljj * Ljj;
          double wj = temp[j];
          double swj2 = beta * wj * wj;
          double gamma = dj * betaPrime + swj2;
          double x = dj + swj2 / betaPrime;
          if (x < 0.0) throw new ArgumentException("Update makes Covariancematrix indefinite");
          double nLjj = Math.Sqrt(x);
          lowerCholesky[j, j] = nLjj;
          betaPrime += swj2 / dj;
          if (j + 1 < n) {
            for (int i = j + 1; i < n; i++) {
              lowerCholesky[i, j] *= a;
            }
            for (int i = j + 1; i < n; i++) {
              temp[i] = wj / Ljj * lowerCholesky[i, j];
            }
            if (gamma == 0) continue;
            for (int i = j + 1; i < n; i++) {
              lowerCholesky[i, j] *= nLjj / Ljj;
            }
            for (int i = j + 1; i < n; i++) {
              lowerCholesky[i, j] += (nLjj * beta * wj / gamma) * temp[i];
            }
          }

        }

      }

      public void mutate(NormalDistributedRandom gauss) {

        //sampling a random z from N(0,I) where I is the Identity matrix;
        int n = lastZ.Length;
        for (int i = 0; i < n; i++) {
          lastZ[i] = gauss.NextDouble();
        }
        //Matrixmultiplication: lastStep = lowerCholesky * lastZ;
        for (int i = 0; i < n; i++) {
          double sum = 0;
          for (int j = 0; j <= i; j++) {
            sum += lowerCholesky[i, j] * lastZ[j];
          }
          lastStep[i] = sum;
        }

        //add the step to x weighted by stepsize;
        for (int i = 0; i < x.Length; i++) {
          x[i] += sigma * lastStep[i];
        }

      }

    }

    //blatantly stolen form HeuristicLab.Optimization.Operators.FastNonDominatedSort
    #region FastNonDominatedSort   
    private enum DominationResult { Dominates, IsDominated, IsNonDominated };

    private List<List<CMAHansenIndividual>> NonDominatedSort(CMAHansenIndividual[] individuals) {
      bool dominateOnEqualQualities = false;
      bool[] maximization = Problem.Maximization;
      if (individuals == null) throw new InvalidOperationException(Name + ": No qualities found.");
      int populationSize = individuals.Length;

      List<List<CMAHansenIndividual>> fronts = new List<List<CMAHansenIndividual>>();
      Dictionary<CMAHansenIndividual, List<int>> dominatedScopes = new Dictionary<CMAHansenIndividual, List<int>>();
      int[] dominationCounter = new int[populationSize];
      //ItemArray<IntValue> rank = new ItemArray<IntValue>(populationSize);

      for (int pI = 0; pI < populationSize - 1; pI++) {
        CMAHansenIndividual p = individuals[pI];
        List<int> dominatedScopesByp;
        if (!dominatedScopes.TryGetValue(p, out dominatedScopesByp))
          dominatedScopes[p] = dominatedScopesByp = new List<int>();
        for (int qI = pI + 1; qI < populationSize; qI++) {
          DominationResult 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) {
            //rank[qI] = new IntValue(0);
            AddToFront(individuals[qI], fronts, 0);
          }
        }
        if (dominationCounter[pI] == 0) {
          //rank[pI] = new IntValue(0);
          AddToFront(p, fronts, 0);
        }
      }
      int i = 0;
      while (i < fronts.Count && fronts[i].Count > 0) {
        List<CMAHansenIndividual> nextFront = new List<CMAHansenIndividual>();
        foreach (CMAHansenIndividual p in fronts[i]) {
          List<int> dominatedScopesByp;
          if (dominatedScopes.TryGetValue(p, out dominatedScopesByp)) {
            for (int k = 0; k < dominatedScopesByp.Count; k++) {
              int dominatedScope = dominatedScopesByp[k];
              dominationCounter[dominatedScope] -= 1;
              if (dominationCounter[dominatedScope] == 0) {
                //rank[dominatedScope] = new IntValue(i + 1);
                nextFront.Add(individuals[dominatedScope]);
              }
            }
          }
        }
        i += 1;
        fronts.Add(nextFront);
      }

      //RankParameter.ActualValue = rank;

      //scope.SubScopes.Clear();

      CMAHansenIndividual[] result = new CMAHansenIndividual[individuals.Length];
      int j = 0;

      for (i = 0; i < fronts.Count; i++) {
        foreach (var p in fronts[i]) {
          p.rank = i;
        }
      }
      return fronts;
    }


    private static void AddToFront(CMAHansenIndividual p, List<List<CMAHansenIndividual>> fronts, int i) {
      if (i == fronts.Count) fronts.Add(new List<CMAHansenIndividual>());
      fronts[i].Add(p);
    }

    private static DominationResult Dominates(CMAHansenIndividual left, CMAHansenIndividual right, bool[] maximizations, bool dominateOnEqualQualities) {
      return Dominates(left.penalizedFitness, right.penalizedFitness, maximizations, dominateOnEqualQualities);
    }

    private static DominationResult Dominates(double[] left, double[] right, 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 (int i = 0; i < left.Length; i++) {
          if (left[i] != right[i]) {
            equal = false;
            break;
          }
        }
        if (equal) return DominationResult.Dominates;
      }

      bool leftIsBetter = false, rightIsBetter = false;
      for (int i = 0; i < left.Length; 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
  }
}