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

using System;
using System.Collections.Generic;
using System.IO;
using System.Linq;
using System.Threading;
using HEAL.Attic;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.Encodings.RealVectorEncoding;
using HeuristicLab.Optimization;
using HeuristicLab.Parameters;
using HeuristicLab.Problems.TestFunctions.MultiObjective;
using HeuristicLab.Problems.TestFunctions.MultiObjective.TestFunctions;
using HeuristicLab.Random;

namespace HeuristicLab.Algorithms.NSGA3
{
    /// <summary>
    /// The Reference Point Based Non-dominated Sorting Genetic Algorithm III was introduced in Deb
    /// et al. 2013. An Evolutionary Many-Objective Optimization Algorithm Using Reference Point
    /// Based Non-dominated Sorting Approach. IEEE Transactions on Evolutionary Computation, 18(4),
    /// pp. 577-601.
    /// </summary>
    [Item("NSGA-III", "The Reference Point Based Non-dominated Sorting Genetic Algorithm III was introduced in Deb et al. 2013. An Evolutionary Many-Objective Optimization Algorithm Using Reference Point Based Non-dominated Sorting Approach. IEEE Transactions on Evolutionary Computation, 18(4), pp. 577-601.")]
    [Creatable(Category = CreatableAttribute.Categories.PopulationBasedAlgorithms, Priority = 136)]
    [StorableType("07C745F7-A8A3-4F99-8B2C-F97E639F9AC3")]
    public class NSGA3 : BasicAlgorithm
    {
        public override bool SupportsPause => true;

        #region ProblemProperties

        public override Type ProblemType
        {
            get { return typeof(MultiObjectiveBasicProblem<RealVectorEncoding>); }
        }

        public new MultiObjectiveBasicProblem<RealVectorEncoding> Problem
        {
            get { return (MultiObjectiveBasicProblem<RealVectorEncoding>)base.Problem; }
            set { base.Problem = value; }
        }

        public int Objectives
        {
            get
            {
                if (!(Problem is MultiObjectiveTestFunctionProblem testFunctionProblem)) throw new NotSupportedException("Only Multi Objective Test Function problems are supported");
                return testFunctionProblem.Objectives;
            }
        }

        #endregion ProblemProperties

        #region Storable fields

        [Storable]
        private IRandom random;

        [Storable]
        private List<Solution> solutions; // maybe todo: rename to nextGeneration (see Run method)

        [Storable]
        private List<ReferencePoint> referencePoints;

        [Storable]
        private NSGA3Selection selection;

        [Storable]
        private double[] allIgds;

        #endregion Storable fields

        #region ParameterAndResultsNames

        // Parameter Names

        private const string PopulationSizeName = "Population Size";
        private const string MaximumGenerationsName = "Maximum Generations";
        private const string CrossoverProbabilityName = "Crossover Probability";
        private const string CrossoverEtaName = "Crossover eta";
        private const string MutationProbabilityName = "Mutation Probability";
        private const string DominateOnEqualQualitiesName = "Dominate On Equal Qualities";
        private const string SetSeedRandomlyName = "Set Seed Randomly";
        private const string SeedName = "Seed";
        private const string AnalyzeEveryGenerationName = "Analyze Every Generation";

        // Results Names

        private const string GeneratedReferencePointsResultName = "Generated Reference Points";
        private const string CurrentGenerationResultName = "Generations";
        private const string GenerationalDistanceResultName = "Generational Distance";
        private const string InvertedGenerationalDistanceResultName = "Inverted Generational Distance";
        private const string HypervolumeResultName = "Hypervolume";
        private const string BestKnownHypervolumeResultName = "Best known hypervolume";
        private const string DifferenceToBestKnownHypervolumeResultName = "Absolute Distance to BestKnownHypervolume";
        private const string ScatterPlotResultName = "Scatter Plot";
        private const string CurrentFrontResultName = "Pareto Front"; // Do not touch this
        private const string RandomSeedUsedName = "Random Seed Used";

        #endregion ParameterAndResultsNames

        #region ParameterProperties

        private IFixedValueParameter<IntValue> PopulationSizeParameter
        {
            get { return (IFixedValueParameter<IntValue>)Parameters[PopulationSizeName]; }
        }

        private IFixedValueParameter<IntValue> MaximumGenerationsParameter
        {
            get { return (IFixedValueParameter<IntValue>)Parameters[MaximumGenerationsName]; }
        }

        private IFixedValueParameter<PercentValue> CrossoverProbabilityParameter
        {
            get { return (IFixedValueParameter<PercentValue>)Parameters[CrossoverProbabilityName]; }
        }

        private IFixedValueParameter<DoubleValue> CrossoverEtaParameter
        {
            get { return (IFixedValueParameter<DoubleValue>)Parameters[CrossoverEtaName]; }
        }

        private IFixedValueParameter<PercentValue> MutationProbabilityParameter
        {
            get { return (IFixedValueParameter<PercentValue>)Parameters[MutationProbabilityName]; }
        }

        private IFixedValueParameter<BoolValue> DominateOnEqualQualitiesParameter
        {
            get { return (IFixedValueParameter<BoolValue>)Parameters[DominateOnEqualQualitiesName]; }
        }

        private IFixedValueParameter<BoolValue> SetSeedRandomlyParameter
        {
            get { return (IFixedValueParameter<BoolValue>)Parameters[SetSeedRandomlyName]; }
        }

        private IFixedValueParameter<IntValue> SeedParameter
        {
            get { return (IFixedValueParameter<IntValue>)Parameters[SeedName]; }
        }

        private IFixedValueParameter<BoolValue> AnalyzeEveryGenerationParameter
        {
            get { return (IFixedValueParameter<BoolValue>)Parameters[AnalyzeEveryGenerationName]; }
        }

        #endregion ParameterProperties

        #region Properties

        public IntValue PopulationSize => PopulationSizeParameter.Value;

        public IntValue MaximumGenerations => MaximumGenerationsParameter.Value;

        public PercentValue CrossoverProbability => CrossoverProbabilityParameter.Value;

        public DoubleValue CrossoverEta => CrossoverEtaParameter.Value;

        public PercentValue MutationProbability => MutationProbabilityParameter.Value;

        public BoolValue DominateOnEqualQualities => DominateOnEqualQualitiesParameter.Value;

        public BoolValue SetSeedRandomly => SetSeedRandomlyParameter.Value;

        public IntValue Seed => SeedParameter.Value;

        public BoolValue AnalyzeEveryGeneration => AnalyzeEveryGenerationParameter.Value;

        #endregion Properties

        #region ResultsProperties

        public DoubleMatrix ResultsGeneratedReferencePoints
        {
            get { return (DoubleMatrix)Results[GeneratedReferencePointsResultName].Value; }
            set { Results[GeneratedReferencePointsResultName].Value = value; }
        }

        public IntValue ResultsCurrentGeneration
        {
            get { return (IntValue)Results[CurrentGenerationResultName].Value; }
            set { Results[CurrentGenerationResultName].Value = value; }
        }

        public DoubleValue ResultsGenerationalDistance
        {
            get { return (DoubleValue)Results[GenerationalDistanceResultName].Value; }
            set { Results[GenerationalDistanceResultName].Value = value; }
        }

        public DoubleValue ResultsInvertedGenerationalDistance
        {
            get { return (DoubleValue)Results[InvertedGenerationalDistanceResultName].Value; }
            set { Results[InvertedGenerationalDistanceResultName].Value = value; }
        }

        public DoubleValue ResultsHypervolume
        {
            get { return (DoubleValue)Results[HypervolumeResultName].Value; }
            set { Results[HypervolumeResultName].Value = value; }
        }

        public DoubleValue ResultsBestKnownHypervolume
        {
            get { return (DoubleValue)Results[BestKnownHypervolumeResultName].Value; }
            set { Results[BestKnownHypervolumeResultName].Value = value; }
        }

        public DoubleValue ResultsDifferenceToBestKnownHypervolume
        {
            get { return (DoubleValue)Results[DifferenceToBestKnownHypervolumeResultName].Value; }
            set { Results[DifferenceToBestKnownHypervolumeResultName].Value = value; }
        }

        public ParetoFrontScatterPlot ResultsScatterPlot
        {
            get { return (ParetoFrontScatterPlot)Results[ScatterPlotResultName].Value; }
            set { Results[ScatterPlotResultName].Value = value; }
        }

        public DoubleMatrix ResultsSolutions
        {
            get { return (DoubleMatrix)Results[CurrentFrontResultName].Value; }
            set { Results[CurrentFrontResultName].Value = value; }
        }

        public IntValue ResultsRandomSeedUsed
        {
            get { return (IntValue)Results[RandomSeedUsedName].Value; }
            set { Results[RandomSeedUsedName].Value = value; }
        }

        #endregion ResultsProperties

        #region Constructors

        public NSGA3() : base()
        {
            Parameters.Add(new FixedValueParameter<IntValue>(PopulationSizeName, "The size of the population of Individuals.", new IntValue(200)));
            Parameters.Add(new FixedValueParameter<IntValue>(MaximumGenerationsName, "The maximum number of generations which should be processed.", new IntValue(1000)));
            Parameters.Add(new FixedValueParameter<PercentValue>(CrossoverProbabilityName, "The probability that the crossover operator is applied on two parents.", new PercentValue(1.0)));
            Parameters.Add(new FixedValueParameter<DoubleValue>(CrossoverEtaName, "TODO: description of this parameter", new DoubleValue(30)));
            Parameters.Add(new FixedValueParameter<PercentValue>(MutationProbabilityName, "The probability that the mutation operator is applied on a Individual.", new PercentValue(0.05)));
            Parameters.Add(new FixedValueParameter<BoolValue>(DominateOnEqualQualitiesName, "Flag which determines wether Individuals with equal quality values should be treated as dominated.", new BoolValue(true)));
            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>(SeedName, "The random seed used to initialize the new pseudo random number generator.", new IntValue(0)));
            Parameters.Add(new FixedValueParameter<BoolValue>(AnalyzeEveryGenerationName, "If set to false, will only calculate the indicators (IGD, etc.) on the last generation and when pausing or stopping the algorithm.", new BoolValue(true)));
        }

        [StorableConstructor]
        public NSGA3(StorableConstructorFlag _) : base(_) { }

        // Each clonable item must have a cloning ctor (deep cloning, the cloner is used to handle
        // cyclic object references). Don't forget to call the cloning ctor of the base class
        public NSGA3(NSGA3 original, Cloner cloner) : base(original, cloner)
        {
            random = cloner.Clone(original.random);
            solutions = original.solutions?.Select(cloner.Clone).ToList();
            referencePoints = original.referencePoints?.Select(cloner.Clone).ToList();
            selection = cloner.Clone(original.selection);
            if (original.allIgds != null)
            {
                allIgds = new double[original.allIgds.Length];
                Array.Copy(original.allIgds, allIgds, allIgds.Length);
            }
        }

        public override IDeepCloneable Clone(Cloner cloner)
        {
            return new NSGA3(this, cloner);
        }

        #endregion Constructors

        #region Initialization

        protected override void Initialize(CancellationToken cancellationToken)
        {
            base.Initialize(cancellationToken);

            InitResults();
            InitFields();
            Analyze(false);
        }

        private void InitResults()
        {
            Results.Add(new Result(GeneratedReferencePointsResultName, "The initially generated reference points", new DoubleMatrix()));
            Results.Add(new Result(CurrentGenerationResultName, "The current generation", new IntValue(1)));
            Results.Add(new Result(GenerationalDistanceResultName, "The generational distance to an optimal pareto front defined in the Problem", new DoubleValue(double.NaN)));
            Results.Add(new Result(InvertedGenerationalDistanceResultName, "The inverted generational distance to an optimal pareto front defined in the Problem", new DoubleValue(double.NaN)));
            Results.Add(new Result(HypervolumeResultName, "The hypervolume of the current front considering the Reference point defined in the Problem", new DoubleValue(double.NaN)));
            Results.Add(new Result(BestKnownHypervolumeResultName, "The best known hypervolume considering the Reference point defined in the Problem", new DoubleValue(double.NaN)));
            Results.Add(new Result(DifferenceToBestKnownHypervolumeResultName, "The difference between the current and the best known hypervolume", new DoubleValue(double.NaN)));
            Results.Add(new Result(ScatterPlotResultName, "A scatterplot displaying the evaluated solutions and (if available) the analytically optimal front", new ParetoFrontScatterPlot()));
            Results.Add(new Result(CurrentFrontResultName, "The Pareto Front", new DoubleMatrix()));
            Results.Add(new Result(RandomSeedUsedName, "The random seed used in the end", new IntValue(-1)));

            if (!(Problem is MultiObjectiveTestFunctionProblem problem)) return;
            var bestKnownFront = problem.BestKnownFront.ToJaggedArray();
            if (problem.BestKnownFront == null) return;
            ResultsScatterPlot = new ParetoFrontScatterPlot(new double[0][], new double[0][], bestKnownFront, problem.Objectives, problem.ProblemSize);
            ResultsBestKnownHypervolume = new DoubleValue(Hypervolume.Calculate(problem.BestKnownFront.ToJaggedArray(), problem.ReferencePoint.CloneAsArray(), problem.Maximization));
        }

        private void InitFields()
        {
            if (SetSeedRandomly.Value)
            {
                System.Random seedRand = new System.Random();
                ResultsRandomSeedUsed.Value = seedRand.Next();
                random = new MersenneTwister(Convert.ToUInt32(ResultsRandomSeedUsed.Value));
            }
            else
            {
                ResultsRandomSeedUsed.Value = Seed.Value;
                random = new MersenneTwister(Convert.ToUInt32(Seed.Value));
            }

            solutions = GetInitialPopulation();
            referencePoints = ReferencePoint.GenerateReferencePoints(random, Objectives);
            ResultsGeneratedReferencePoints = Utility.ConvertToDoubleMatrix(referencePoints);
            var problem = Problem as MultiObjectiveTestFunctionProblem;

            selection = new NSGA3Selection(problem.Objectives, problem.Maximization, PopulationSize.Value, random, DominateOnEqualQualities.Value);

            allIgds = new double[MaximumGenerations.Value];
        }

        private List<Solution> GetInitialPopulation()
        {
            var problem = Problem as MultiObjectiveTestFunctionProblem;
            if (problem.Bounds.Rows != 1) throw new Exception();

            // Initialise solutions
            List<Solution> solutions = new List<Solution>(PopulationSize.Value);
            for (int i = 0; i < PopulationSize.Value; i++)
            {
                double minBound = problem.Bounds[0, 0];
                double maxBound = problem.Bounds[0, 1];
                RealVector randomRealVector = new RealVector(Problem.Encoding.Length, random, minBound, maxBound);
                var solution = new Solution(randomRealVector);
                solution.Objectives = Evaluate(solution.Chromosome);
                solutions.Add(solution);
            }

            return solutions;
        }

        #endregion Initialization

        #region Overriden Methods

        protected override void Run(CancellationToken cancellationToken)
        {
            try
            {
                while (ResultsCurrentGeneration.Value < MaximumGenerations.Value)
                {
                    List<Solution> qt = Mutate(Recombine(solutions));
                    foreach (var solution in qt)
                    {
                        solution.Objectives = Evaluate(solution.Chromosome);
                    }

                    List<Solution> rt = Utility.Concat(solutions, qt);

                    solutions = selection.SelectSolutionsForNextGeneration(rt, GetCopyOfReferencePoints());

                    Analyze(!AnalyzeEveryGeneration.Value);

                    ResultsCurrentGeneration.Value++;
                    cancellationToken.ThrowIfCancellationRequested();
                }
            }
            catch (OperationCanceledException ex)
            {
                throw new OperationCanceledException("Optimization process was cancelled.", ex);
            }
            catch (Exception ex)
            {
                throw new Exception($"Failed in generation {ResultsCurrentGeneration}.", ex);
            }
            finally
            {
                Analyze(false);
                WriteIGDsToFile();
            }
        }

        private void WriteIGDsToFile()
        {
            string resultsDirectory = "Results";
            Directory.CreateDirectory(resultsDirectory);
            var problem = Problem as MultiObjectiveTestFunctionProblem;
            string fileName = Path.Combine(resultsDirectory, $"{problem.TestFunction.Name}({problem.Objectives})_DOEQ={DominateOnEqualQualities.Value}_Seed={ResultsRandomSeedUsed.Value}.txt");

            File.WriteAllLines(fileName, allIgds.Select(igd => igd.ToString()));
        }

        #endregion Overriden Methods

        #region Private Methods

        private List<ReferencePoint> GetCopyOfReferencePoints()
        {
            if (referencePoints == null) return null;

            return referencePoints.Select(rp => new ReferencePoint(rp)).ToList();
        }

        private void Analyze(bool calculateOnlyIGD)
        {
            if (!calculateOnlyIGD)
            {
                ResultsScatterPlot = new ParetoFrontScatterPlot(solutions.Select(x => x.Objectives).ToArray(), solutions.Select(x => x.Chromosome.ToArray()).ToArray(), ResultsScatterPlot.ParetoFront, ResultsScatterPlot.Objectives, ResultsScatterPlot.ProblemSize);
                ResultsSolutions = solutions.Select(s => s.Chromosome.ToArray()).ToMatrix();
            }

            if (!(Problem is MultiObjectiveTestFunctionProblem problem)) return;

            CalculateInvertedGenerationalDistance(problem);

            if (!calculateOnlyIGD)
            {
                ResultsGenerationalDistance = new DoubleValue(problem.BestKnownFront != null ? GenerationalDistance.Calculate(solutions.Select(s => s.Objectives), problem.BestKnownFront.ToJaggedArray(), 1) : double.NaN);
                var front = NonDominatedSelect.GetDominatingVectors(solutions.Select(x => x.Objectives), problem.ReferencePoint.CloneAsArray(), Problem.Maximization, true).ToArray();
                if (front.Length == 0) return;
                ResultsHypervolume = new DoubleValue(Hypervolume.Calculate(front, problem.ReferencePoint.CloneAsArray(), problem.Maximization));
                ResultsDifferenceToBestKnownHypervolume = new DoubleValue(ResultsBestKnownHypervolume.Value - ResultsHypervolume.Value);

                Problem.Analyze(
                    solutions.Select(s => (Individual)new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, s.Chromosome) } })).ToArray(),
                    solutions.Select(s => s.Objectives).ToArray(),
                    Results,
                    random);
            }
        }

        private void CalculateInvertedGenerationalDistance(MultiObjectiveTestFunctionProblem problem)
        {
            double igd;
            if (problem.TestFunction is ScaledDTLZ dtlz)
            {
                // Normalize objectives before calculating the igd
                var scalingFactor = dtlz.GetScalingFactor(problem.Objectives);
                var front = new double[solutions.Count][];

                for (int i = 0; i < solutions.Count; i++)
                {
                    front[i] = new double[problem.Objectives];
                    Array.Copy(solutions[i].Objectives, front[i], front[i].Length);
                    for (int objective = 0; objective < problem.Objectives; objective++)
                    {
                        // get normalized value
                        front[i][objective] = front[i][objective] / Math.Pow(scalingFactor, objective);
                    }
                }

                // Normalize best known front before calculating the igd
                var bestKnownFront = problem.BestKnownFront.ToJaggedArray();
                foreach (var point in bestKnownFront)
                {
                    for (int objective = 0; objective < problem.Objectives; objective++)
                    {
                        point[objective] = point[objective] / Math.Pow(scalingFactor, objective);
                    }
                }

                igd = InvertedGenerationalDistance.Calculate(front, bestKnownFront, 1);
            }
            else
            {
                igd = problem.BestKnownFront != null ? InvertedGenerationalDistance.Calculate(solutions.Select(s => s.Objectives), problem.BestKnownFront.ToJaggedArray(), 1) : double.NaN;
            }
            ResultsInvertedGenerationalDistance = new DoubleValue(igd);
            allIgds[ResultsCurrentGeneration.Value - 1] = ResultsInvertedGenerationalDistance.Value;
        }

        /// <summary>
        /// Returns the fitness of the given <paramref name="chromosome" /> by applying the Evaluate
        /// method of the Problem.
        /// </summary>
        /// <param name="chromosome"></param>
        /// <returns></returns>
        private double[] Evaluate(RealVector chromosome)
        {
            return Problem.Evaluate(new SingleEncodingIndividual(Problem.Encoding, new Scope { Variables = { new Variable(Problem.Encoding.Name, chromosome) } }), random);
        }

        private List<Solution> Recombine(List<Solution> solutions)
        {
            List<Solution> childSolutions = new List<Solution>();

            for (int i = 0; i < solutions.Count; i += 2)
            {
                if (solutions.Count != PopulationSize.Value) throw new Exception("solutions.Count != PopulationSize");
                int parentIndex1 = random.Next(solutions.Count);
                int parentIndex2 = random.Next(solutions.Count);
                // ensure that the parents are not the same object
                if (parentIndex1 == parentIndex2) parentIndex2 = (parentIndex2 + 1) % solutions.Count;
                var parent1 = solutions[parentIndex1];
                var parent2 = solutions[parentIndex2];

                // Do crossover with crossoverProbabilty == 1 in order to guarantee that a crossover happens
                var children = SimulatedBinaryCrossover.Apply(random,
                    Problem.Encoding.Bounds, parent1.Chromosome, parent2.Chromosome, CrossoverProbability.Value, CrossoverEta.Value);

                childSolutions.Add(new Solution(children.Item1));
                childSolutions.Add(new Solution(children.Item2));
            }

            return childSolutions;
        }

        private List<Solution> Mutate(List<Solution> solutions, double eta = 20)
        {
            foreach (var solution in solutions)
            {
                for (int i = 0; i < solution.Chromosome.Length; i++)
                {
                    if (random.NextDouble() > MutationProbability.Value) continue;

                    double y = solution.Chromosome[i];
                    double lb;
                    double ub;
                    var problem = Problem as MultiObjectiveTestFunctionProblem;
                    if (problem.Bounds.Rows == 1) lb = problem.Bounds[0, 0];
                    else lb = problem.Bounds[i, 0];
                    if (problem.Bounds.Rows == 1) ub = problem.Bounds[0, 1];
                    else ub = problem.Bounds[i, 1];

                    double delta1 = (y - lb) / (ub - lb);
                    double delta2 = (ub - y) / (ub - lb);

                    double mut_pow = 1.0 / (eta + 1.0);

                    double rnd = random.NextDouble();
                    double deltaq;
                    if (rnd <= 0.5)
                    {
                        double xy = 1.0 - delta1;
                        double val = 2.0 * rnd + (1.0 - 2.0 * rnd) * (Math.Pow(xy, (eta + 1.0)));
                        deltaq = Math.Pow(val, mut_pow) - 1.0;
                    }
                    else
                    {
                        double xy = 1.0 - delta2;
                        double val = 2.0 * (1.0 - rnd) + 2.0 * (rnd - 0.5) * (Math.Pow(xy, (eta + 1.0)));
                        deltaq = 1.0 - (Math.Pow(val, mut_pow));
                    }

                    y += deltaq * (ub - lb);
                    y = Math.Min(ub, Math.Max(lb, y));

                    solution.Chromosome[i] = y;
                }
            }
            return solutions;
        }

        #endregion Private Methods
    }
}