source: branches/3057_DynamicALPS/HeuristicLab.Algorithms.DynamicALPS/3.4/DynamicALPSAlgorithmBase.cs @ 17729

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2019 Heuristic and Evolutionary Algorithms Laboratory (HEAL)
 *
 * 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

// 03/02/2020
// TODO LIST:                                                                     
// 1. Dynamic reference point strategy                                               
// 2. Normalized fitness value strategy, desibility function.                         
// 3. HVC calculation should be definitely improved, at least in the 2D case.
// 4. multiple point strategy when $\lambda>1$ 

using HEAL.Attic;
using HeuristicLab.Analysis;
using HeuristicLab.Common;
using HeuristicLab.Core;
using HeuristicLab.Data;
using HeuristicLab.ExpressionGenerator;
using HeuristicLab.Optimization;
using HeuristicLab.Parameters;
using HeuristicLab.Problems.DataAnalysis;
using HeuristicLab.Problems.TestFunctions.MultiObjective;
using HeuristicLab.Random;
using System;
using System.Collections.Generic;
using System.Drawing;
using System.Linq;
using CancellationToken = System.Threading.CancellationToken;

namespace HeuristicLab.Algorithms.DynamicALPS {
  [Item("DynamicALPSAlgorithmBase", "Base class for all DynamicALPS algorithm variants.")]
  [StorableType("C0141748-DF5A-4CA0-A3DD-1DFB98236F7E")]
  public abstract class DynamicALPSAlgorithmBase : BasicAlgorithm {
    #region data members
 
    [StorableType("75C9EA99-D699-4A1F-8AB2-AB86B7A2FD68")]
    protected enum NeighborType { NEIGHBOR, POPULATION }


    [StorableType("2A71E397-05CE-460F-B982-EE2F4B37C354")]
    // TCHE = Chebyshev (Tchebyshev)
    // PBI  = Penalty-based boundary intersection
    // AGG  = Weighted sum
    public enum FunctionType { TCHE, PBI, AGG }

    [Storable]
    protected double[] IdealPoint { get; set; }
    [Storable]
    protected double[] NadirPoint { get; set; } // potentially useful for objective normalization

    [Storable]
    protected double[][] lambda_moead;

    [Storable]
    protected int[][] neighbourhood;

    [Storable]
    protected IDynamicALPSSolution[] solutions;

    [Storable]
    protected FunctionType functionType;

    [Storable]
    protected IDynamicALPSSolution[] population;

    [Storable]
    protected IDynamicALPSSolution[][] layerPopulation;

    [Storable]
    protected bool[] activeLayer;

    [Storable]
    protected double[] layerCrossoverProbability;

    [Storable]
    protected IDynamicALPSSolution[][] layerDiscardPopulation;

    [Storable]
    protected IDynamicALPSSolution[] layerDiscardIndivdual;
    

    [Storable]
    protected IDynamicALPSSolution[][] layerOffspringPopulation;

    [Storable]
    protected IDynamicALPSSolution[][] layerJointPopulation;

    [Storable]
    protected IDynamicALPSSolution[] offspringPopulation;

    //[Storable]
    //protected IDynamicALPSSolution[] jointPopulation;

    [Storable]
    protected int evaluatedSolutions;

    [Storable]
    protected ExecutionContext executionContext;

    [Storable]
    protected IScope globalScope;

    [Storable]
    protected ExecutionState previousExecutionState;

    [Storable]
    protected ExecutionState executionState;

    protected DoubleArray ReferencePoint {
      get {
        if (Problem is MultiObjectiveTestFunctionProblem) {
          var problem = (MultiObjectiveTestFunctionProblem)Problem;
          return problem.ReferencePoint;
        }
        else {
          return null;
        }
      }
    }
    #endregion

    #region parameters
    private const string SeedParameterName = "Seed";
    private const string SetSeedRandomlyParameterName = "SetSeedRandomly";
    private const string PopulationSizeParameterName = "PopulationSize";
    private const string ResultPopulationSizeParameterName = "ResultPopulationSize";
    private const string CrossoverProbabilityParameterName = "CrossoverProbability";
    private const string CrossoverParameterName = "Crossover";
    private const string MutationProbabilityParameterName = "MutationProbability";
    private const string MutatorParameterName = "Mutator";
    private const string MaximumEvaluatedSolutionsParameterName = "MaximumEvaluatedSolutions";
    private const string RandomParameterName = "Random";
    private const string AnalyzerParameterName = "Analyzer";
    
    
    // MOEA-D parameters
    //private const string NeighbourSizeParameterName = "NeighbourSize";
    //private const string NeighbourhoodSelectionProbabilityParameterName = "NeighbourhoodSelectionProbability";
    //private const string MaximumNumberOfReplacedSolutionsParameterName = "MaximumNumberOfReplacedSolutions";
    //private const string FunctionTypeParameterName = "FunctionType";
    // private const string NormalizeObjectivesParameterName = "NormalizeObjectives";

    // SMS-EMOA parameters: 
    private const string LambdaParameterName = "Lambda";   // The number of offspring size 
    private const string ALPSLayersParameterName = "ALPSLayers";   // The number of offspring size 
    private const string ALPSAgeGapParameterName = "ALPSAgeGap";   // The number of offspring size 
    private const string InitializeLayerPopulationMethodName = "InitializationLayerPopulations";



    // "Parameters" are defined in "HeuristicLab.Parameters"
    // Contains: generic parameters of every class/algorithm/instance, 
    // It seems that "I***ValueParameter" is declared in "Heuristic.core", where "***ValueParameter" are defined in "HeuristicLab.Parameter"
    // The function of "I***ValueParameter" is to bridge current scripts to "HeuristicLab.Parameter". 
    public IValueParameter<MultiAnalyzer> AnalyzerParameter {
      get { return (ValueParameter<MultiAnalyzer>)Parameters[AnalyzerParameterName]; }
    }

    //public IConstrainedValueParameter<StringValue> FunctionTypeParameter
    //{
    //  get { return (IConstrainedValueParameter<StringValue>)Parameters[FunctionTypeParameterName]; }
    //}
    //public IFixedValueParameter<IntValue> NeighbourSizeParameter
    //{
    //  get { return (IFixedValueParameter<IntValue>)Parameters[NeighbourSizeParameterName]; }
    //}
    //public IFixedValueParameter<BoolValue> NormalizeObjectivesParameter
    //{
    //  get { return (IFixedValueParameter<BoolValue>)Parameters[NormalizeObjectivesParameterName]; }
    //}
    //public IFixedValueParameter<IntValue> MaximumNumberOfReplacedSolutionsParameter
    //{
    //  get { return (IFixedValueParameter<IntValue>)Parameters[MaximumNumberOfReplacedSolutionsParameterName]; }
    //}
    //public IFixedValueParameter<DoubleValue> NeighbourhoodSelectionProbabilityParameter
    //{
    //  get { return (IFixedValueParameter<DoubleValue>)Parameters[NeighbourhoodSelectionProbabilityParameterName]; }
    //}
    public IFixedValueParameter<IntValue> SeedParameter {
      get { return (IFixedValueParameter<IntValue>)Parameters[SeedParameterName]; }
    }
    public IFixedValueParameter<BoolValue> SetSeedRandomlyParameter {
      get { return (IFixedValueParameter<BoolValue>)Parameters[SetSeedRandomlyParameterName]; }
    }
    private IValueParameter<IntValue> PopulationSizeParameter {
      get { return (IValueParameter<IntValue>)Parameters[PopulationSizeParameterName]; }
    }
    // KF, SMS-EMOA
    private IValueParameter<IntValue> LambdaParameter {
      get { return (IValueParameter<IntValue>)Parameters[LambdaParameterName]; }
    }
    //// KF, DynamicALPS
    private IValueParameter<IntValue> ALPSLayersParameter{
      get { return (IValueParameter<IntValue>)Parameters[ALPSLayersParameterName]; }
    }
    private IValueParameter<IntValue> ALPSAgeGapParameter {
      get { return (IValueParameter<IntValue>)Parameters[ALPSAgeGapParameterName]; }
    }

    private IValueParameter<BoolValue> ALPSInitialzeLayerPopulationParameter {
      get { return (IValueParameter<BoolValue>)Parameters[InitializeLayerPopulationMethodName]; }
    }
    

    private IValueParameter<IntValue> ResultPopulationSizeParameter {
      get { return (IValueParameter<IntValue>)Parameters[ResultPopulationSizeParameterName]; }
    }

    public IValueParameter<PercentValue> CrossoverProbabilityParameter {
      get { return (IValueParameter<PercentValue>)Parameters[CrossoverProbabilityParameterName]; }
    }
    public IConstrainedValueParameter<ICrossover> CrossoverParameter {
      get { return (IConstrainedValueParameter<ICrossover>)Parameters[CrossoverParameterName]; }
    }
    public IValueParameter<PercentValue> MutationProbabilityParameter {
      get { return (IValueParameter<PercentValue>)Parameters[MutationProbabilityParameterName]; }
    }
    public IConstrainedValueParameter<IManipulator> MutatorParameter {
      get { return (IConstrainedValueParameter<IManipulator>)Parameters[MutatorParameterName]; }
    }
    public IValueParameter<IntValue> MaximumEvaluatedSolutionsParameter {
      get { return (IValueParameter<IntValue>)Parameters[MaximumEvaluatedSolutionsParameterName]; }
    }
    public IValueParameter<IRandom> RandomParameter {
      get { return (IValueParameter<IRandom>)Parameters[RandomParameterName]; }
    }
    #endregion

    #region parameter properties
    public new IMultiObjectiveHeuristicOptimizationProblem Problem {
      get { return (IMultiObjectiveHeuristicOptimizationProblem)base.Problem; }
      set { base.Problem = 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 IntValue PopulationSize {
      get { return PopulationSizeParameter.Value; }
      set { PopulationSizeParameter.Value = value; }
    }
    public IntValue Lambda {
      get { return LambdaParameter.Value; }
      set { LambdaParameter.Value = value; }
    }

    public IntValue ResultPopulationSize {
      get { return ResultPopulationSizeParameter.Value; }
      set { ResultPopulationSizeParameter.Value = value; }
    }

    public IntValue ALPSLayers {
      get { return ALPSLayersParameter.Value; }
      set { ALPSLayersParameter.Value = value; }
    }

    public IntValue ALPSAgeGap {
      get { return ALPSAgeGapParameter.Value; }
      set { ALPSAgeGapParameter.Value = value; }
    }
    public BoolValue UseAverageAge {
      get { return ALPSInitialzeLayerPopulationParameter.Value; }
      set { ALPSInitialzeLayerPopulationParameter.Value = value; }
    }

    public PercentValue CrossoverProbability {
      get { return CrossoverProbabilityParameter.Value; }
      set { CrossoverProbabilityParameter.Value = value; }
    }
    public ICrossover Crossover {
      get { return CrossoverParameter.Value; }
      set { CrossoverParameter.Value = value; }
    }
    public PercentValue MutationProbability {
      get { return MutationProbabilityParameter.Value; }
      set { MutationProbabilityParameter.Value = value; }
    }
    public IManipulator Mutator {
      get { return MutatorParameter.Value; }
      set { MutatorParameter.Value = value; }
    }
    public MultiAnalyzer Analyzer {
      get { return AnalyzerParameter.Value; }
      set { AnalyzerParameter.Value = value; }
    }
    public IntValue MaximumEvaluatedSolutions {
      get { return MaximumEvaluatedSolutionsParameter.Value; }
      set { MaximumEvaluatedSolutionsParameter.Value = value; }
    }
    #endregion

    #region constructors
    public DynamicALPSAlgorithmBase() {
      // Add or define or specify the parameters that may be use in SMS-EMOA.   
      // ***("Name", "Description", "Value")
      //  Name                            Type                Description
      //  FixedValueParameter:            ANY                 Not changed???
      //  ValueParameter:                                     Changable??? What is the difference between "ValueParameter" and "FixedVlaueParameter"?????


      // types: 
      //      IntValue
      //      BoolValue
      //      DoubleValue
      //      PercentValue
      //      ICrossover:       
      //      IManipulator:     
      //      IRandom:         
      //      MultiAnalyzer:    
      //      ---------
      Parameters.Add(new FixedValueParameter<IntValue>(SeedParameterName, "The random seed used to initialize the new pseudo random number generator.", new IntValue(0)));
      Parameters.Add(new FixedValueParameter<BoolValue>(SetSeedRandomlyParameterName, "True if the random seed should be set to a random value, otherwise false.", new BoolValue(true)));
      Parameters.Add(new ValueParameter<IntValue>(PopulationSizeParameterName, "The size of the population of solutions.", new IntValue(100)));
      Parameters.Add(new ValueParameter<IntValue>(ResultPopulationSizeParameterName, "The size of the population of solutions.", new IntValue(100)));
      Parameters.Add(new ValueParameter<PercentValue>(CrossoverProbabilityParameterName, "The probability that the crossover operator is applied.", new PercentValue(0.9)));
      Parameters.Add(new ConstrainedValueParameter<ICrossover>(CrossoverParameterName, "The operator used to cross solutions."));
      Parameters.Add(new ValueParameter<PercentValue>(MutationProbabilityParameterName, "The probability that the mutation operator is applied on a solution.", new PercentValue(0.25)));
      Parameters.Add(new ConstrainedValueParameter<IManipulator>(MutatorParameterName, "The operator used to mutate solutions."));
      Parameters.Add(new ValueParameter<MultiAnalyzer>("Analyzer", "The operator used to analyze each generation.", new MultiAnalyzer()));
      Parameters.Add(new ValueParameter<IntValue>(MaximumEvaluatedSolutionsParameterName, "The maximum number of evaluated solutions (approximately).", new IntValue(100_000)));
      Parameters.Add(new ValueParameter<IRandom>(RandomParameterName, new FastRandom()));
      Parameters.Add(new ValueParameter<BoolValue>(InitializeLayerPopulationMethodName, "Whether use average age to initialize the layer population or not. If not, move the older individuals to layer populations", new BoolValue(true)));

      // SMS-EMOA, kf
      Parameters.Add(new ValueParameter<IntValue>(LambdaParameterName, "The size of the offsprings. Now, it only works when lambda = 1", new IntValue(1)));
      // DynamicALPS, KF
      Parameters.Add(new ValueParameter<IntValue>(ALPSLayersParameterName, "Test, maximum = 1000, defualt is 1", new IntValue(10)));
      Parameters.Add(new ValueParameter<IntValue>(ALPSAgeGapParameterName, "Test, maximum = 1000, defualt is 20", new IntValue(20)));
      

    }

    protected DynamicALPSAlgorithmBase(DynamicALPSAlgorithmBase original, Cloner cloner) : base(original, cloner) {
      functionType = original.functionType;
      evaluatedSolutions = original.evaluatedSolutions;
      previousExecutionState = original.previousExecutionState;

      if (original.IdealPoint != null) {
        IdealPoint = (double[])original.IdealPoint.Clone();
      }

      if (original.NadirPoint != null) {
        NadirPoint = (double[])original.NadirPoint.Clone();
      }

      if (original.lambda_moead != null) {
        lambda_moead = (double[][])original.lambda_moead.Clone();
      }

      if (original.neighbourhood != null) {
        neighbourhood = (int[][])original.neighbourhood.Clone();
      }

      if (original.solutions != null) {
        solutions = original.solutions.Select(cloner.Clone).ToArray();
      }

      if (original.population != null) {
        population = original.population.Select(cloner.Clone).ToArray();
      }

      if (original.offspringPopulation != null) {
        offspringPopulation = original.offspringPopulation.Select(cloner.Clone).ToArray();
      }

      //if (original.jointPopulation != null) {
      //  jointPopulation = original.jointPopulation.Select(x => cloner.Clone(x)).ToArray();
      //}

      if (original.executionContext != null) {
        executionContext = cloner.Clone(original.executionContext);
      }

      if (original.globalScope != null) {
        globalScope = cloner.Clone(original.globalScope);
      }
    }



    [StorableConstructor]
    protected DynamicALPSAlgorithmBase(StorableConstructorFlag deserializing) : base(deserializing) { }
    #endregion

    private void InitializePopulation(ExecutionContext executionContext, CancellationToken cancellationToken, IRandom random, bool[] maximization) {
      // creator: how to create the initilized population. "UniformRandom" is used here. 
      // TODO: LHS, latin hypercube sampling? Exisit??? 
      var creator = Problem.SolutionCreator;
      var evaluator = Problem.Evaluator;

      // dimensions: objective space
      var dimensions = maximization.Length;
      var populationSize = PopulationSize.Value;
      population = new IDynamicALPSSolution[populationSize];

      var parentScope = executionContext.Scope;
      // first, create all individuals
      for (int i = 0; i < populationSize; ++i) {
        var childScope = new Scope(i.ToString()) { Parent = parentScope };
        ExecuteOperation(executionContext, cancellationToken, executionContext.CreateChildOperation(creator, childScope));
        parentScope.SubScopes.Add(childScope);
      }

      for (int i = 0; i < populationSize; ++i) {
        var childScope = parentScope.SubScopes[i];
        ExecuteOperation(executionContext, cancellationToken, executionContext.CreateChildOperation(evaluator, childScope));

        var qualities = (DoubleArray)childScope.Variables["Qualities"].Value;

        // solution: a method, contains a decision vector and objecitve values     
        //    solution.Qualities:     objective values, fitness values 
        //    solution.Individual:    decision vector
        var solution = new DynamicALPSSolution(childScope, dimensions, 0);
        for (int j = 0; j < dimensions; ++j) {
          // TODO: convert maximization problems into minimization problems. 
          solution.Qualities[j] = maximization[j] ? 1 - qualities[j] : qualities[j];
        }

        // population is a collection of solution.  
        population[i] = solution;

        // KF, DyanmicALPS
        population[i].HypervolumeContribution[0] = -0;
        population[i].NondominanceRanking[0] = -0;
        population[i].Age = 1;
        population[i].IndividualPc = CrossoverProbability.Value;
        population[i].IndividualPc = MutationProbability.Value;
      }
    }

    protected void InitializeAlgorithm(CancellationToken cancellationToken) { // Type of random operator, "FastRandom" in this script. 
      // RandomParameter <-- Parameters in "HeuristicLab.Core.ParameterizedNameItem", 
      var rand = RandomParameter.Value;

      // Initialize random seed
      // If random seed exist, get it; otherwise, 
      if (SetSeedRandomly) Seed = RandomSeedGenerator.GetSeed();

      // Call 
      rand.Reset(Seed);

      bool[] maximization = ((BoolArray)Problem.MaximizationParameter.ActualValue).CloneAsArray();

      // dimensions: the dimension in an objective space 
      var dimensions = maximization.Length;


      var populationSize = PopulationSize.Value;

      InitializePopulation(executionContext, cancellationToken, rand, maximization);

      IdealPoint = new double[dimensions];
      IdealPoint.UpdateIdeal(population);

      NadirPoint = Enumerable.Repeat(double.MinValue, dimensions).ToArray();
      //NadirPoint = new double[dimensions];
      NadirPoint.UpdateNadir(population);


      evaluatedSolutions = populationSize;
    }

    protected override void Initialize(CancellationToken cancellationToken) {
      globalScope = new Scope("Global Scope");
      executionContext = new ExecutionContext(null, this, globalScope);

      // set the execution context for parameters to allow lookup
      foreach (var parameter in Problem.Parameters.OfType<IValueParameter>()) {
        // we need all of these in order for the wiring of the operators to work
        globalScope.Variables.Add(new Variable(parameter.Name, parameter.Value));
      }
      globalScope.Variables.Add(new Variable("Results", Results)); // make results available as a parameter for analyzers etc.

      base.Initialize(cancellationToken);
    }

    public override bool SupportsPause => true;




    // Mate Selection. 
    // Randomly select a specific number of individuals for later operators. 
    // Inputs: 
    //    1. random:                      Random number generate method
    //    2. numberOfSolutionToSelect:    The number of selection    
    // Outputs: 
    //    1. listOfSolutions:             The selection individuals
    protected List<int> MatingSelection(IRandom random, int numberOfSolutionsToSelect) {
      int populationSize = PopulationSize.Value;

      var listOfSolutions = new List<int>(numberOfSolutionsToSelect);

      while (listOfSolutions.Count < numberOfSolutionsToSelect) {
        var selectedSolution = random.Next(populationSize);

        bool flag = true;
        foreach (int individualId in listOfSolutions) {
          if (individualId == selectedSolution) {
            flag = false;
            break;
          }
        }

        if (flag) {
          listOfSolutions.Add(selectedSolution);
        }
      }
      return listOfSolutions;
    }

    protected void ApplyCrossover(int lambda) {
    }

    protected void ApplyMutation(int lambda) {
    }


    protected void ApplyEvaluation(int lambda) {
    }

    protected void ApplyMateSelection(int rho) {
    }

    protected void InitializeLayer(int indexLayer, int populationSize, int lambda) {
      layerPopulation[indexLayer] = new IDynamicALPSSolution[populationSize];
      layerJointPopulation[indexLayer] = new IDynamicALPSSolution[populationSize + lambda];
      layerOffspringPopulation[indexLayer] = new IDynamicALPSSolution[lambda];
      layerDiscardPopulation[indexLayer] = new IDynamicALPSSolution[populationSize];
      activeLayer[indexLayer] = true;
    }


    // Select/Discard the individual(s) according to HVC 
    protected void SmetricSelection(int lambda, int nLayerALPS) {
      var wholePopulation = layerJointPopulation[nLayerALPS];
      var qualities = wholePopulation.Select(x => x.Qualities).ToArray();

      var maxPoint = Enumerable.Range(0, IdealPoint.Length).Select(idx => qualities.Max(q => q[idx])).ToArray();

      var maximization = Enumerable.Repeat(false, IdealPoint.Length).ToArray(); // Minimization or maximization ????
      var pf2 = DominationCalculator<IDynamicALPSSolution>.CalculateAllParetoFronts(wholePopulation, qualities, maximization, out int[] ranking);

      int numberOfLayer;             // number of layers in PF
      int numberOfLastLayer;          // number of discarded points in PF (the number of points in the last layer)

      pf2.RemoveAt(pf2.Count() - 1);
      numberOfLayer = pf2.Count();
      numberOfLastLayer = pf2[numberOfLayer - 1].Count();
      double[] hvc = new double[numberOfLastLayer];
      int discardIndex;
      if (numberOfLastLayer > lambda) {
        double tempHV;
        double smetric;
        var lastLayer = pf2.Last();

        // TODO: This can be use for dynamic reference point strategy later. Kaifeng , 02/2020
        // smetric = Hypervolume.Calculate(lastLayer.Select(x => x.Item2), Enumerable.Repeat(11d, NadirPoint.Length).ToArray(), maximization);

        var reference = Enumerable.Repeat(double.MaxValue, maximization.Length).ToArray(); // TODO Dynamic Reference point for each layer. Maximum * 1.1
                                                                                           //if (nLayerALPS != 0) {
        for (int i = 0; i < reference.Length; i++) {
          reference[i] = 1.1 * maxPoint[i];
          if (reference[i] > 10000) {
            reference[i] = 9999;          // set a upper bound for the reference point
          }
        }
        //}
        //else {
        //  reference = ReferencePoint.ToArray();
        //}

        var nondominated = NonDominatedSelect.GetDominatingVectors(lastLayer.Select(x => x.Item2), reference, maximization, false);
        smetric = nondominated.Any() ? Hypervolume.Calculate(nondominated, reference, maximization) : int.MinValue;

        for (int ii = 0; ii < lastLayer.Count; ++ii) {
          try { // TODO: This can be use for dynamic reference point strategy later. Kaifeng , 02/2020
            // tempHV = Hypervolume.Calculate(indices.Where(idx => idx != ii).Select(idx => lastLayer[idx].Item2), Enumerable.Repeat(11d, NadirPoint.Length).ToArray(), maximization);
            tempHV = Hypervolume.Calculate(Enumerable.Range(0, lastLayer.Count).Where(idx => idx != ii).Select(idx => lastLayer[idx].Item2), reference, maximization);
          }
          catch {
            tempHV = int.MinValue;
          }
          hvc[ii] = smetric - tempHV;
          tempHV = 0;
        }

        discardIndex = Array.IndexOf(hvc, hvc.Min());
        //layerDiscardPopulation[nLayerALPS] = pf2[numberOfLayer - 1][discardIndex].Item1.ToArray();
        layerDiscardIndivdual[nLayerALPS] = pf2[numberOfLayer - 1].Select(x => x.Item1).ToArray()[discardIndex];
        pf2[numberOfLayer - 1].RemoveAt(discardIndex);
      }
      else {
        // TODO: This should be updated when $lambda > 1$
        discardIndex = pf2.Count() - 1;
        layerDiscardIndivdual[nLayerALPS] = pf2[discardIndex].Select(x => x.Item1).ToArray()[0];
        pf2.RemoveAt(pf2.Count() - 1);
        numberOfLayer = numberOfLayer - 1;
        
      }
      layerPopulation[nLayerALPS] = pf2.SelectMany(x => x.Select(y => y.Item1)).ToArray();
      
    }

    public static double SampleGaussian(IRandom random, double mean, double stddev) {
      // The method requires sampling from a uniform random of (0,1]
      // but Random.NextDouble() returns a sample of [0,1).
      double x1 = 1 - random.NextDouble();
      double x2 = 1 - random.NextDouble();

      double y1 = Math.Sqrt(-2.0 * Math.Log(x1)) * Math.Cos(2.0 * Math.PI * x2);
      return y1 * stddev + mean;
    }

    protected int SMSEMOA(int populationSize, int lambda, int counterLayerALPS) {
      var innerToken = new CancellationToken();
      bool[] maximization = ((BoolArray)Problem.MaximizationParameter.ActualValue).CloneAsArray();
      var maximumEvaluatedSolutions = MaximumEvaluatedSolutions.Value;
      var crossover = Crossover;
      var crossoverProbability = layerCrossoverProbability[0];
      var mutator = Mutator;
      var mutationProbability = MutationProbability.Value;
      var evaluator = Problem.Evaluator;
      var analyzer = Analyzer;
      var rand = RandomParameter.Value;


      int indexOffspring = 0;
      var mates = MatingSelection(rand, 2); // select parents
                                            //var s1 = (IScope)population[mates[0]].Individual.Clone();
                                            //var s2 = (IScope)population[mates[1]].Individual.Clone();
                                            //var ages = population.Select(x => x.Age).ToArray();

      var s1 = (IScope)layerPopulation[counterLayerALPS][mates[0]].Individual.Clone();
      var s2 = (IScope)layerPopulation[counterLayerALPS][mates[1]].Individual.Clone();
      var ages = layerPopulation[counterLayerALPS].Select(x => x.Age).ToArray();

      var s1_age = ages[mates[0]];
      var s2_age = ages[mates[1]];
      int offSpringAge = 0;
      s1.Parent = s2.Parent = globalScope;
      IScope childScope = null;

      // crossoverProbability = crossoverProbability - 0.02* counterLayerALPS;
      //var test = SampleGaussian(rand, 0, 1);

      //crossoverProbability = 1 / (1 + Math.Exp(-0.02 * SampleGaussian(rand, 0, 1)) * (1 - crossoverProbability) / crossoverProbability);

      if (crossoverProbability < 0.5)
        crossoverProbability = 0.5;
      else if(crossoverProbability > 0.95)
      {
        crossoverProbability = 0.95;
      }
      layerCrossoverProbability[counterLayerALPS] = crossoverProbability;

      // crossover
      if (rand.NextDouble() < crossoverProbability) {
        childScope = new Scope($"{mates[0]}+{mates[1]}") { Parent = executionContext.Scope };
        childScope.SubScopes.Add(s1);
        childScope.SubScopes.Add(s2);
        var opCrossover = executionContext.CreateChildOperation(crossover, childScope);
        ExecuteOperation(executionContext, innerToken, opCrossover);
        offSpringAge = Math.Max(s1_age, s2_age) + 1;
        childScope.SubScopes.Clear(); // <<-- VERY IMPORTANT!
      }
      else { // MUTATION   POLISHI
        if (childScope == null) {
          offSpringAge = ages[mates[0]];
        }
        else {
        }
        childScope = childScope ?? s1;
        var opMutation = executionContext.CreateChildOperation(mutator, childScope);
        ExecuteOperation(executionContext, innerToken, opMutation);
        offSpringAge = offSpringAge + 1;
      }
      if (childScope != null) { // Evaluate the childScope
        var opEvaluation = executionContext.CreateChildOperation(evaluator, childScope);
        ExecuteOperation(executionContext, innerToken, opEvaluation);
        // childScope 
        var qualities = (DoubleArray)childScope.Variables["Qualities"].Value;
        var childSolution = new DynamicALPSSolution(childScope, maximization.Length, 0);
        // set child qualities
        for (int j = 0; j < maximization.Length; ++j) {
          childSolution.Qualities[j] = maximization[j] ? 1 - qualities[j] : qualities[j];
        }
        IdealPoint.UpdateIdeal(childSolution.Qualities);
        NadirPoint.UpdateNadir(childSolution.Qualities);
        // TODO, KF -- For later usage when $lambda > 1$
        childSolution.HypervolumeContribution = null;
        childSolution.NondominanceRanking = null;
        childSolution.Age = offSpringAge;
        layerOffspringPopulation[counterLayerALPS][indexOffspring] = childSolution;
        ++evaluatedSolutions;
        indexOffspring += 1;
      }
      else {
        // no crossover or mutation were applied, a child was not produced, do nothing
      }


      layerJointPopulation[counterLayerALPS] = new IDynamicALPSSolution[populationSize + lambda];
      layerPopulation[counterLayerALPS].CopyTo(layerJointPopulation[counterLayerALPS], 0);
      layerOffspringPopulation[counterLayerALPS].CopyTo(layerJointPopulation[counterLayerALPS], populationSize);

      SmetricSelection(lambda, counterLayerALPS);   // SMS-EMOA 
      return evaluatedSolutions;
    }






      // Update the Pareto-front approximation set and scatter the solutions in PF approximation set. 
      protected ResultCollection UpdateParetoFronts(IDynamicALPSSolution[] solutions, double[] IdealPoint) {
 

      //var qualities = population.Select(x => Enumerable.Range(0, NadirPoint.Length).Select(i => x.Qualities[i] / NadirPoint[i]).ToArray()).ToArray();
      var qualities = solutions.Select(x => x.Qualities).ToArray();
      var maximization = Enumerable.Repeat(false, IdealPoint.Length).ToArray();                             // DynamicALPS minimizes everything internally
      var pf = DominationCalculator<IDynamicALPSSolution>.CalculateBestParetoFront(solutions, qualities, maximization);

      var pf2 = DominationCalculator<IDynamicALPSSolution>.CalculateAllParetoFronts(solutions, qualities, maximization, out int[] ranking);
      var n = (int)EnumerableExtensions.BinomialCoefficient(IdealPoint.Length, 2);


      // Struture hypervolume 
      // [0,0]:  Value of HV
      // [0,1]:  PF size, $|PF|$
      var hypervolumes = new DoubleMatrix(n == 1 ? 1 : n + 1, 2) { ColumnNames = new[] { "PF hypervolume", "PF size" } };


      // HV calculation 
      // pf.Select(x => x.Item2): the "Item2" in var "pd"
      // Enumerable.Repeat(1d, NadirPoint.Length).ToArray():     reference point 
      // maximization:   type of optimization problem:
      //               True:  maximization problem 
      //               False: minimization problem
      var reference = Enumerable.Repeat(double.MaxValue, maximization.Length).ToArray();
      if (ReferencePoint is null) {     // KF, 20200217 -- fix no reference point on real-world applications. If No reference points in Algorithms, use 1.1 \times max objective values as the reference point
        for (int i = 0; i < reference.Length; i++) {
          reference[i] = 1.1 * reference[i];
          if (reference[i] > 10000) {
            reference[i] = 9999;          // set a upper bound for the reference point
          }
        }
      }
      else {
        reference = ReferencePoint.ToArray();
      }

      var nondominated = NonDominatedSelect.GetDominatingVectors(pf.Select(x => x.Item2), reference, maximization, false);
      hypervolumes[0, 0] = nondominated.Any() ? Hypervolume.Calculate(nondominated, reference, maximization) : int.MinValue;

      //hypervolumes[0, 0] = Hypervolume.Calculate(pf.Select(x => x.Item2), reference, maximization);
      hypervolumes[0, 1] = pf.Count;
      Console.WriteLine("Current HV is", hypervolumes[0, 0]);

      var elementNames = new List<string>() { "Pareto Front" };
      var results = new ResultCollection();

      ResultCollection innerResults;
      if (results.ContainsKey("Hypervolume Analysis")) {
        innerResults = (ResultCollection)results["Hypervolume Analysis"].Value;
      }
      else {
        innerResults = new ResultCollection();
        results.AddOrUpdateResult("Hypervolume Analysis", innerResults);
      }

      ScatterPlot sp;
      if (IdealPoint.Length == 2) {
        var points = pf.Select(x => new Point2D<double>(x.Item2[0], x.Item2[1]));
        var r = OnlinePearsonsRCalculator.Calculate(points.Select(x => x.X), points.Select(x => x.Y), out OnlineCalculatorError error);
        if (error != OnlineCalculatorError.None) { r = double.NaN; }
        var resultName = "Pareto Front Analysis ";
        if (!innerResults.ContainsKey(resultName)) {
          sp = new ScatterPlot() {
            //VisualProperties = {
            //  XAxisMinimumAuto = true, XAxisMinimumFixedValue = 0d, XAxisMaximumAuto = false, XAxisMaximumFixedValue = 1d,
            //  YAxisMinimumAuto = true, YAxisMinimumFixedValue = 0d, YAxisMaximumAuto = false, YAxisMaximumFixedValue = 1d
            //}
          };
          sp.Rows.Add(new ScatterPlotDataRow(resultName, "", points) { VisualProperties = { PointSize = 8 } });
          innerResults.AddOrUpdateResult(resultName, sp);
        }
        else {
          sp = (ScatterPlot)innerResults[resultName].Value;
          sp.Rows[resultName].Points.Replace(points);
        }
        sp.Name = $"Dimensions [0, 1], correlation: {r.ToString("N2")}";
      }
      else if (IdealPoint.Length > 2) {
        var indices = Enumerable.Range(0, IdealPoint.Length).ToArray();
        var visualProperties = new ScatterPlotDataRowVisualProperties { PointSize = 8, Color = Color.LightGray };
        var combinations = indices.Combinations(2).ToArray();
        var maximization2d = new[] { false, false };
        var solutions2d = pf.Select(x => x.Item1).ToArray();
        for (int i = 0; i < combinations.Length; ++i) {
          var c = combinations[i].ToArray();

          // calculate the hypervolume in the 2d coordinate space
          var reference2d = new[] { 1d, 1d };
          var qualities2d = pf.Select(x => new[] { x.Item2[c[0]], x.Item2[c[1]] }).ToArray();
          var pf2d = DominationCalculator<IDynamicALPSSolution>.CalculateBestParetoFront(solutions2d, qualities2d, maximization2d);

          hypervolumes[i + 1, 0] = pf2d.Count > 0 ? Hypervolume.Calculate(pf2d.Select(x => x.Item2), reference2d, maximization2d) : 0d;
          hypervolumes[i + 1, 1] = pf2d.Count;

          var resultName = $"Pareto Front Analysis [{c[0]}, {c[1]}]";
          elementNames.Add(resultName);

          var points = pf.Select(x => new Point2D<double>(x.Item2[c[0]], x.Item2[c[1]]));
          var pf2dPoints = pf2d.Select(x => new Point2D<double>(x.Item2[0], x.Item2[1]));

          if (!innerResults.ContainsKey(resultName)) {
            sp = new ScatterPlot() {
              VisualProperties = {
                XAxisMinimumAuto = false, XAxisMinimumFixedValue = 0d, XAxisMaximumAuto = false, XAxisMaximumFixedValue = 1d,
                YAxisMinimumAuto = false, YAxisMinimumFixedValue = 0d, YAxisMaximumAuto = false, YAxisMaximumFixedValue = 1d
              }
            };
            sp.Rows.Add(new ScatterPlotDataRow("Pareto Front", "", points) { VisualProperties = visualProperties });
            sp.Rows.Add(new ScatterPlotDataRow($"Pareto Front [{c[0]}, {c[1]}]", "", pf2dPoints) { VisualProperties = { PointSize = 10, Color = Color.OrangeRed } });
            innerResults.AddOrUpdateResult(resultName, sp);
          }
          else {
            sp = (ScatterPlot)innerResults[resultName].Value;
            sp.Rows["Pareto Front"].Points.Replace(points);
            sp.Rows[$"Pareto Front [{c[0]}, {c[1]}]"].Points.Replace(pf2dPoints);
          }
          var r = OnlinePearsonsRCalculator.Calculate(points.Select(x => x.X), points.Select(x => x.Y), out OnlineCalculatorError error);
          var r2 = r * r;
          sp.Name = $"Pareto Front [{c[0]}, {c[1]}], correlation: {r2.ToString("N2")}";
        }
      }
      hypervolumes.RowNames = elementNames;
      innerResults.AddOrUpdateResult("Hypervolumes", hypervolumes);

      return results;
    }

    #region operator wiring and events
    protected void ExecuteOperation(ExecutionContext executionContext, CancellationToken cancellationToken, IOperation operation) {
      Stack<IOperation> executionStack = new Stack<IOperation>();
      executionStack.Push(operation);
      while (executionStack.Count > 0) {
        cancellationToken.ThrowIfCancellationRequested();
        IOperation next = executionStack.Pop();
        if (next is OperationCollection) {
          OperationCollection coll = (OperationCollection)next;
          for (int i = coll.Count - 1; i >= 0; i--)
            if (coll[i] != null) executionStack.Push(coll[i]);
        }
        else if (next is IAtomicOperation) {
          IAtomicOperation op = (IAtomicOperation)next;
          next = op.Operator.Execute((IExecutionContext)op, cancellationToken);
          if (next != null) executionStack.Push(next);
        }
      }
    }

    protected virtual void UpdateAnalyzers() {
      Analyzer.Operators.Clear();
      if (Problem != null) {
        foreach (IAnalyzer analyzer in Problem.Operators.OfType<IAnalyzer>()) {
          foreach (IScopeTreeLookupParameter param in analyzer.Parameters.OfType<IScopeTreeLookupParameter>())
            param.Depth = 1;
          Analyzer.Operators.Add(analyzer, analyzer.EnabledByDefault);
        }
      }
    }

    private void UpdateCrossovers() {
      ICrossover oldCrossover = CrossoverParameter.Value;
      CrossoverParameter.ValidValues.Clear();
      ICrossover defaultCrossover = Problem.Operators.OfType<ICrossover>().FirstOrDefault();

      foreach (ICrossover crossover in Problem.Operators.OfType<ICrossover>().OrderBy(x => x.Name))
        CrossoverParameter.ValidValues.Add(crossover);

      if (oldCrossover != null) {
        ICrossover crossover = CrossoverParameter.ValidValues.FirstOrDefault(x => x.GetType() == oldCrossover.GetType());
        if (crossover != null) CrossoverParameter.Value = crossover;
        else oldCrossover = null;
      }
      if (oldCrossover == null && defaultCrossover != null)
        CrossoverParameter.Value = defaultCrossover;
    }

    private void UpdateMutators() {
      IManipulator oldMutator = MutatorParameter.Value;
      MutatorParameter.ValidValues.Clear();
      IManipulator defaultMutator = Problem.Operators.OfType<IManipulator>().FirstOrDefault();

      foreach (IManipulator mutator in Problem.Operators.OfType<IManipulator>().OrderBy(x => x.Name))
        MutatorParameter.ValidValues.Add(mutator);

      if (oldMutator != null) {
        IManipulator mutator = MutatorParameter.ValidValues.FirstOrDefault(x => x.GetType() == oldMutator.GetType());
        if (mutator != null) MutatorParameter.Value = mutator;
        else oldMutator = null;
      }

      if (oldMutator == null && defaultMutator != null)
        MutatorParameter.Value = defaultMutator;
    }

    protected override void OnProblemChanged() {
      UpdateCrossovers();
      UpdateMutators();
      UpdateAnalyzers();
      base.OnProblemChanged();
    }

    protected override void OnExecutionStateChanged() {
      previousExecutionState = executionState;
      executionState = ExecutionState;
      base.OnExecutionStateChanged();
    }

    public void ClearState() {
      solutions = null;
      population = null;
      offspringPopulation = null;
      //jointPopulation = null;
      lambda_moead = null;
      neighbourhood = null;
      if (executionContext != null && executionContext.Scope != null) {
        executionContext.Scope.SubScopes.Clear();
      }
    }

    protected override void OnStopped() {
      ClearState();
      base.OnStopped();
    }
    #endregion
  }
}