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

using System;
using System.Collections.Generic;
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.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;

        #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 MutationProbabilityName = "Mutation Probability";
        private const string DominateOnEqualQualitiesName = "Dominate On Equal Qualities";
        private const string SetSeedRandomlyName = "Set Seed Randomly";
        private const string SeedName = "Seed";

        // 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

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

        #endregion ParameterProperties

        #region Properties

        public IntValue PopulationSize => PopulationSizeParameter.Value;

        public IntValue MaximumGenerations => MaximumGenerationsParameter.Value;

        public PercentValue CrossoverProbability => CrossoverProbabilityParameter.Value;

        public PercentValue MutationProbability => MutationProbabilityParameter.Value;

        public BoolValue DominateOnEqualQualities => DominateOnEqualQualitiesParameter.Value;

        public BoolValue SetSeedRandomly => SetSeedRandomlyParameter.Value;
        public IntValue Seed => SeedParameter.Value;

        // todo: create one property for the Generated Reference Points and one for the current
        // generations reference points

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

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

        // Persistence uses this ctor to improve deserialization efficiency. If we would use the
        // default ctor instead this would completely initialize the object (e.g. creating
        // parameters) even though the data is later overwritten by the stored data.
        [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();
        }

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

        #endregion Constructors

        #region Initialization

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

            SetParameters();
            InitResults();
            InitFields();
            Analyze();
        }

        // todo: trigger this when the problem parameters change
        private void SetParameters()
        {
            // Set population size
            int numberOfGeneratedReferencePoints = ReferencePoint.GetNumberOfGeneratedReferencePoints(Objectives);
            if (numberOfGeneratedReferencePoints == -1) throw new NotSupportedException("The number of objectives is not supported.");
            PopulationSize.Value = ReferencePoint.GetPopulationSizeForReferencePoints(numberOfGeneratedReferencePoints);

            // Set mutation probability
            MutationProbability.Value = 1.0 / Problem.Encoding.Length;

            // todo: Set MaximumGenerations in Problem definition
        }

        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()));

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

        private void InitFields()
        {
            random = new MersenneTwister();
            if (!SetSeedRandomly.Value)
                random.Reset(Seed.Value);

            solutions = GetInitialPopulation();
            referencePoints = ReferencePoint.GenerateReferencePoints(random, Objectives);
            ResultsGeneratedReferencePoints = Utility.ConvertToDoubleMatrix(referencePoints);
        }

        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.Fitness = 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)
                {
                    // todo: make parameter out of this
                    List<Solution> qt = Mutate(Recombine(solutions), 30);
                    foreach (var solution in qt)
                        solution.Fitness = Evaluate(solution.Chromosome);

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

                    // todo: remove this check
                    for (int i = 0; i < rt.Count / 2; i++)
                        if (!solutions.Contains(rt[i])) throw new Exception($"This should never happen: !solutions.Contains(rt[{i}])");

                    solutions = NSGA3Selection.SelectSolutionsForNextGeneration(rt, GetCopyOfReferencePoints(), Problem.Maximization, PopulationSize.Value, random);

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

        #endregion Overriden Methods

        #region Private Methods

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

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

        private void Analyze()
        {
            ResultsScatterPlot = new ParetoFrontScatterPlot(solutions.Select(x => x.Fitness).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;

            // Indicators
            ResultsGenerationalDistance = new DoubleValue(problem.BestKnownFront != null ? GenerationalDistance.Calculate(solutions.Select(s => s.Fitness), problem.BestKnownFront.ToJaggedArray(), 1) : double.NaN);
            ResultsInvertedGenerationalDistance = new DoubleValue(problem.BestKnownFront != null ? InvertedGenerationalDistance.Calculate(solutions.Select(s => s.Fitness), problem.BestKnownFront.ToJaggedArray(), 1) : double.NaN);

            var front = NonDominatedSelect.GetDominatingVectors(solutions.Select(x => x.Fitness), 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);
        }

        /// <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)
            {
                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);

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

            return childSolutions;
        }

        private List<Solution> Mutate(List<Solution> solutions, double eta)
        {
            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
    }
}