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

#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.Random;
using HeuristicLab.Problems.MultiObjectiveTestFunctions;
using System.Linq;

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]; }
        }
        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; }
        }

        #endregion

        #region ResultsProperties
        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; }
        }

        //Datatable
        private DataTable ResultsQualities
        {
            get { return ((DataTable)Results["Timetable"].Value); }
        }
        private DataRow ResultsHypervolumeBest
        {
            get { return ResultsQualities.Rows["Best Hypervolume"]; }
        }
        private DataRow ResultsHypervolumeIteration
        {
            get { return ResultsQualities.Rows["Iteration Hypervolume"]; }
        }
        private DataRow ResultsGenerationalDistanceIteration
        {
            get { return ResultsQualities.Rows["Iteration Generational Distance"]; }
        }
        private DataRow ResultsInvertedGenerationalDistanceIteration
        {
            get { return ResultsQualities.Rows["Iteration Inverted Generational Distance"]; }
        }
        private DataRow ResultsCrowdingIteration
        {
            get { return ResultsQualities.Rows["Iteration Crowding"]; }
        }


        //QualityIndicators
        private double ResultsHypervolume
        {
            get { return ((DoubleValue)Results["Hypervolume"].Value).Value; }
            set { ((DoubleValue)Results["Hypervolume"].Value).Value = value; }
        }
        private double ResultsGenerationalDistance
        {
            get { return ((DoubleValue)Results["Generational Distance"].Value).Value; }
            set { ((DoubleValue)Results["Generational Distance"].Value).Value = value; }
        }
        private double ResultsInvertedGenerationalDistance
        {
            get { return ((DoubleValue)Results["Inverted Generational Distance"].Value).Value; }
            set { ((DoubleValue)Results["Inverted Generational Distance"].Value).Value = value; }
        }
        private double ResultsCrowding
        {
            get { return ((DoubleValue)Results["Crowding"].Value).Value; }
            set { ((DoubleValue)Results["Crowding"].Value).Value = value; }
        }
        private double ResultsBestHypervolume
        {
            get { return ((DoubleValue)Results["Best Hypervolume"].Value).Value; }
            set { ((DoubleValue)Results["Best Hypervolume"].Value).Value = value; }
        }
        private double ResultsBestKnownHypervolume
        {
            get { return ((DoubleValue)Results["Best Known Hypervolume"].Value).Value; }
            set { ((DoubleValue)Results["Best Known Hypervolume"].Value).Value = value; }
        }

        //Solutions
        private DoubleMatrix ResultsSolutions
        {
            get { return ((DoubleMatrix)Results["Solutions"].Value); }
            set { Results["Solutions"].Value = value; }
        }
        private MOSolution ResultsSolutionsPlot
        {
            get { return ((MOSolution)Results["Solutions Scatterplot"].Value); }
            set { Results["Solutions Scatterplot"].Value = value; }
        }
        #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 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, Problem);   //TODO: This is a battel in its own right to be fought another day
                front_[lc].selected = false;
                front_.Swap(lc, front_.Count - 1);
                front_.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);
                ResultsEvaluations++;
                offspring.penalizedFitness = offspring.fitness;
            }
            else
            {
                RealVector t = closestFeasible(offspring.x);
                offspring.fitness = Problem.Evaluate(t, random);
                ResultsEvaluations++;
                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++)
            {
                int 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)
        {
            DoubleMatrix bounds = Problem.Bounds;
            for (int i = 0; i < offspring.Length; i++)
            {
                int dim = i % bounds.Rows;
                if (bounds[dim, 0] > offspring[i] || offspring[i] > bounds[dim, 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 < PopulationSize; 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()
        {
            solutions = new CMAHansenIndividual[PopulationSize];
            for (int i = 0; i < PopulationSize; i++)
            {
                RealVector x = new RealVector(Problem.ProblemSize); // uniform distibution in all dimesions assumed TODO is there a better way to initialize RandomVectors
                var bounds = Problem.Bounds;
                for (int j = 0; j < Problem.Objectives; j++)
                {
                    int 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()
        {
            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;

        }


        private void initResults()
        {
            Results.Add(new Result("Iterations", new IntValue(0)));
            Results.Add(new Result("Evaluations", new IntValue(0)));

            Results.Add(new Result("Hypervolume", new DoubleValue(0)));
            Results.Add(new Result("Best Hypervolume", new DoubleValue(0)));
            Results.Add(new Result("Generational Distance", new DoubleValue(0)));
            Results.Add(new Result("Inverted Generational Distance", new DoubleValue(0)));
            Results.Add(new Result("Crowding", new DoubleValue(0)));

            var table = new DataTable("QualityIndicators");
            Results.Add(new Result("Timetable", table));
            table.Rows.Add(new DataRow("Best Hypervolume"));
            table.Rows.Add(new DataRow("Iteration Hypervolume"));
            table.Rows.Add(new DataRow("Iteration Crowding"));
            table.Rows.Add(new DataRow("Iteration Generational Distance"));
            table.Rows.Add(new DataRow("Iteration Inverted Generational Distance"));



            Results.Add(new Result("Solutions", new DoubleMatrix()));
            Results.Add(new Result("Solutions Scatterplot", new MOSolution(null, null, Problem.BestKnownFront.ToArray<double[]>(), Problem.Objectives)));


            if (Problem.BestKnownFront != null)
            {
                Results.Add(new Result("Best Known Hypervolume", new DoubleValue(Hypervolume.Calculate(Problem.BestKnownFront, Problem.TestFunction.ReferencePoint(Problem.Objectives), Problem.Maximization))));
            }
        }
        #endregion

        #region analyze
        private void analyzeSolutions()
        {
            //solutions
            ResultsSolutionsPlot = new MOSolution(solutions.Select(x => x.fitness).ToArray<double[]>(),
                solutions.Select(x => ToArray(x.x)).ToArray<double[]>(),
                ResultsSolutionsPlot.ParetoFront,
                ResultsSolutionsPlot.Objectives);
            ResultsSolutions = ToMatrix(solutions.Select(x => ToArray(x.x)));
            analyzeQualityIndicators();

        }

        private void analyzeQualityIndicators()
        {

            var front = NonDominatedSelect.selectNonDominatedVectors(solutions.Select(x => x.fitness), Problem.Maximization, true);
            front = NonDominatedSelect.removeNonReferenceDominatingVectors(front, Problem.TestFunction.ReferencePoint(Problem.Objectives), Problem.Maximization, false);
            var bounds = ToArray(Problem.Bounds);
            ResultsCrowding = Crowding.Calculate(front, bounds);
            ResultsGenerationalDistance = Problem.BestKnownFront != null ? GenerationalDistance.Calculate(front, Problem.BestKnownFront, 1) : Double.NaN;

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

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

            //Datalines
            ResultsHypervolumeBest.Values.Add(ResultsBestHypervolume);
            ResultsHypervolumeIteration.Values.Add(ResultsHypervolume);
            ResultsCrowdingIteration.Values.Add(ResultsCrowding);
            ResultsGenerationalDistanceIteration.Values.Add(ResultsGenerationalDistance);
            ResultsInvertedGenerationalDistanceIteration.Values.Add(ResultsInvertedGenerationalDistance);


        }
        #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);

            initResults();
            initStrategy();
            initSolutions();

            // Loop until iteration limit reached or canceled.
            for (ResultsIterations = 1; ResultsIterations < MaximumGenerations; ResultsIterations++)
            {
                try
                {
                    iterate();
                    cancellationToken.ThrowIfCancellationRequested();
                }
                finally
                {
                    analyzeSolutions();
                }
            }
        }

        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);
            var parents = merge(solutions, offspring);
            selection(parents, solutions.Length);
            updatePopulation(parents);

        }
        #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.evolutionPath = (RealVector)other.evolutionPath.Clone();
                this.x = (RealVector)other.x.Clone();


                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;
                lastZ = new RealVector(x.Length);
                int n = lastZ.Length;
                for (int i = 0; i < n; i++)
                {
                    lastZ[i] = gauss.NextDouble();
                }
                //Matrixmultiplication: lastStep = lowerCholesky * lastZ;
                lastStep = new RealVector(x.Length);
                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];

            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



        #region conversions

        public double[] ToArray(RealVector r)
        {
            double[] d = new double[r.Length];
            for (int i = 0; i < r.Length; i++)
            {
                d[i] = r[i];
            }
            return d;
        }

        public DoubleMatrix ToMatrix(IEnumerable<double[]> data)
        {
            var d2 = data.ToArray<double[]>();
            DoubleMatrix mat = new DoubleMatrix(d2.Length, d2[0].Length);
            for (int i = 0; i < mat.Rows; i++)
            {
                for (int j = 0; j < mat.Columns; j++)
                {
                    mat[i, j] = d2[i][j];
                }
            }
            return mat;
        }
        public double[,] ToArray(DoubleMatrix data)
        {

            double[,] mat = new double[data.Rows, data.Columns];
            for (int i = 0; i < data.Rows; i++)
            {
                for (int j = 0; j < data.Columns; j++)
                {
                    mat[i, j] = data[i, j];
                }
            }
            return mat;
        }

        #endregion
    }
}