Changeset 17479 for branches/3057_DynamicALPS/HeuristicLab.Algorithms.DynamicALPS

Timestamp:

03/17/20 16:29:48 (5 years ago)

Author:

kyang

Message:

#3057

1. upload the latest version of ALPS with SMS-EMOA
2. upload the related dynamic test problems (dynamic, single-objective symbolic regression), written by David Daninel.

Location:

branches/3057_DynamicALPS/HeuristicLab.Algorithms.DynamicALPS/3.4

Files:

• DynamicALPSASolution.cs (modified) (1 diff)
• DynamicALPSAlgorithm.cs (modified) (10 diffs)
• DynamicALPSAlgorithmBase.cs (modified) (21 diffs)

branches/3057_DynamicALPS/HeuristicLab.Algorithms.DynamicALPS/3.4/DynamicALPSASolution.cs

@@ -115 +115 @@
     {
       Qualities = (double[])original.Qualities.Clone();
-      Constraints = (double[])original.Qualities.Clone();
-      HypervolumeContribution = (double[])original.HypervolumeContribution.Clone();       // kf, sms-emoa
-      NondominanceRanking = (int[])original.NondominanceRanking.Clone();                  // kf, sms-emoa
-      Individual = (IItem)original.Individual.Clone(cloner);
+      Constraints = (double[])original.Constraints.Clone();
+      if (original.HypervolumeContribution != null) {
+        HypervolumeContribution = (double[])original.HypervolumeContribution.Clone();       // kf, sms-emoa
+      }
+      if (original.NondominanceRanking != null) {
+        NondominanceRanking = (int[])original.NondominanceRanking.Clone();                  // kf, sms-emoa
+      }
+      if (original.Individual != null) {
+        Individual = (IItem)original.Individual.Clone(cloner);
+      }
     }
 

branches/3057_DynamicALPS/HeuristicLab.Algorithms.DynamicALPS/3.4/DynamicALPSAlgorithm.cs

@@ -120 +120 @@
 
       // IMPROVE: ........
-      int[] ageGapArray = new int[] { 10, 20, 30, 40, 50, 60, 70, 80, 90, 100 };
+      //                              1   2   4   8   16  32  64  128 256 512
+      int[] ageGapArray = new int[] { 20, 40, 80, 160, 320, 640, 1280, 2560, 5120, 10240 };
       int[] numberDiscard = new int[10];
 
@@ -126 +127 @@
 
       activeLayer = new bool[nLayerALPS];
+      layerCrossoverProbability = new double[nLayerALPS];
       int[][] ageMatrix = new int[nLayerALPS][];      // layer * population size
 
@@ -132 +134 @@
       // cancellation token for the inner operations which should not be immediately cancelled
       var innerToken = new CancellationToken();
-
-
-      // run some analyzer on each layer (for now calculate scatter plots )
-      List<Tuple<int, ScatterPlotAnalyzer>> layerAnalyzers;
 
 
@@ -143 +141 @@
       layerJointPopulation = new IDynamicALPSSolution[nLayerALPS][];
       layerDiscardPopulation = new IDynamicALPSSolution[nLayerALPS][];
+      layerDiscardIndivdual = new IDynamicALPSSolution[nLayerALPS];
 
 
@@ -148 +147 @@
       // BUG: The size of offspring should vary in different layers!!!!
       layerOffspringPopulation[counterLayerALPS] = new IDynamicALPSSolution[lambda];
+      // layerDiscardPopulation[counterLayerALPS] = new IDynamicALPSSolution[populationSize];  // for the previous version, is used to store the individuals whose age is older than the age gap
       layerDiscardPopulation[counterLayerALPS] = new IDynamicALPSSolution[populationSize];
       population.CopyTo(layerPopulation[counterLayerALPS], 0);
 
       activeLayer[counterLayerALPS] = true;
+      layerCrossoverProbability[counterLayerALPS] = 0;
+      var test = UseAverageAge.Value;
       for (int i = 0; i < nLayerALPS; i++) {
         if (i == 0) {
@@ -158 +160 @@
         else { activeLayer[i] = false; }
         numberDiscard[i] = 0;
+        layerCrossoverProbability[i] = CrossoverProbability.Value;
       }
-
+      int bottomLayerIDX = 0;
+      int godScope = 0;
       // Mainloop 
       while (evaluatedSolutions < maximumEvaluatedSolutions && !cancellationToken.IsCancellationRequested) {
         for (int i = 0; i < nLayerALPS; i++) {            // loop for every layer
           int discardedIndividualIndex = 0;
-          if (activeLayer[i] == true) {                   // check the layer is active or not. 
-            evaluatedSolutions = SMSEMOA(populationSize, lambda, i);
+          var currentLayerIDX = i;
+          var nextLayerIDX = i + 1;
+          if (nextLayerIDX == nLayerALPS) {
+            nextLayerIDX = bottomLayerIDX;
+            godScope = 1;
+          }
+          else { godScope = 0; }
+
+
+          if (activeLayer[currentLayerIDX] == true) {                   // check the current layer is active or not. 
+            evaluatedSolutions = SMSEMOA(populationSize, lambda, currentLayerIDX);                  // get the offspring -- layerJointPopulation
             if (evaluatedSolutions >= maximumEvaluatedSolutions) { break; }           // check evaluation budget
-            ageMatrix[i] = layerJointPopulation[i].Select(x => x.Age).ToArray();      // get age info of the current layer joint population 
-            if (ageMatrix[i].Max() > ageGapArray[i]) {    // mature: moving 
-              discardedIndividualIndex = ageMatrix[i].ToList().IndexOf(ageMatrix[i].Max());   // BUG when two individual has the same maximal age???? NOT POSSBILE IN SMS-EMOA 
-              layerDiscardPopulation[i][numberDiscard[i]] = layerJointPopulation[i][discardedIndividualIndex];                    // move the individual to the next layer 
-              layerJointPopulation[i].Where((x, idx) => idx != discardedIndividualIndex).ToArray().CopyTo(layerPopulation[i], 0); // discard the indivudal in the current layer  
-              numberDiscard[i] += 1;
-              if (activeLayer[i + 1] == false) {   // next layer is not active  // bug, if i == number of layer,  out of range .
-                if (numberDiscard[i] == populationSize) {
-                  InitializeLayer(i + 1, populationSize, lambda);  // initilizae the layerPopulation for the next layer && ACTIVE FLAGE
-                  layerDiscardPopulation[i].CopyTo(layerPopulation[i + 1], 0);
-                  numberDiscard[i] = 0;
+            ageMatrix[currentLayerIDX] = layerJointPopulation[currentLayerIDX].Select(x => x.Age).ToArray();      // get age info of the current layer joint population 
+            #region version 1: use average to initialize the layer population
+            if (UseAverageAge.Value == true) {
+              if (activeLayer[nextLayerIDX] == false) {// next layer is not active yet
+                if (ageMatrix[currentLayerIDX].Average() > ageGapArray[currentLayerIDX]) {    // the next layer is initialized 
+                  InitializeLayer(nextLayerIDX, populationSize, lambda);  // initilizae the layerPopulation for the next layer && ACTIVE FLAGE
+                  layerJointPopulation[currentLayerIDX].CopyTo(layerJointPopulation[nextLayerIDX], 0);
+                  SmetricSelection(lambda, nextLayerIDX);                //   layerpopulation is updated here, 
+                }
+                else {// the next layer is not initialized 
+                  SmetricSelection(lambda, currentLayerIDX);
+                }
+              }
+              else { // next layer is active 
+                if (activeLayer.All(x => x) && godScope == 1) {  // all the layers are active and the current layer is the top layer, move the discarded individual from the top to bottom, and reset the age
+                  SmetricSelection(lambda, currentLayerIDX);
+                  layerPopulation[bottomLayerIDX].CopyTo(layerJointPopulation[bottomLayerIDX], 0);
+                  layerDiscardIndivdual[currentLayerIDX].Age = 0;
+                  layerJointPopulation[bottomLayerIDX][populationSize] = layerDiscardIndivdual[currentLayerIDX];
+                  SmetricSelection(lambda, bottomLayerIDX);
+                }
+                else {
+                  if (ageMatrix[currentLayerIDX].Max() > ageGapArray[currentLayerIDX]) {      // moving 
+                    discardedIndividualIndex = ageMatrix[currentLayerIDX].ToList().IndexOf(ageMatrix[currentLayerIDX].Max());
+                    layerPopulation[nextLayerIDX].CopyTo(layerJointPopulation[nextLayerIDX], 0);
+                    layerJointPopulation[currentLayerIDX].Where((x, idx) => idx == discardedIndividualIndex).ToArray().CopyTo(layerJointPopulation[nextLayerIDX], populationSize);
+                    SmetricSelection(lambda, nextLayerIDX); // AGE and HVC 
+                    layerJointPopulation[currentLayerIDX].Where((x, idx) => idx != discardedIndividualIndex).ToArray().CopyTo(layerPopulation[currentLayerIDX], 0);  // dicard the individual in the current layer
+                  }
+                  else {   // next layer is active, but the age is not mature.
+                    SmetricSelection(lambda, currentLayerIDX);
+                  }
+                }
+              }
+              #endregion
+            }
+            else { 
+            #region version 2: use individual age to to initialize the next layer
+            if (ageMatrix[currentLayerIDX].Max() > ageGapArray[currentLayerIDX]) {    // mature: moving 
+              discardedIndividualIndex = ageMatrix[currentLayerIDX].ToList().IndexOf(ageMatrix[currentLayerIDX].Max());   // BUG when two individual has the same maximal age???? NOT POSSBILE IN SMS-EMOA 
+              layerDiscardPopulation[currentLayerIDX][numberDiscard[currentLayerIDX]] = layerJointPopulation[currentLayerIDX][discardedIndividualIndex];                    // move the individual to the next layer 
+              layerJointPopulation[currentLayerIDX].Where((x, idx) => idx != discardedIndividualIndex).ToArray().CopyTo(layerPopulation[currentLayerIDX], 0); // discard the indivudal in the current layer  
+              numberDiscard[currentLayerIDX] += 1;
+              if (activeLayer[nextLayerIDX] == false) {   // next layer is not active  // bug, if i == number of layer,  out of range .
+                if (numberDiscard[currentLayerIDX] == populationSize) {
+                  InitializeLayer(nextLayerIDX, populationSize, lambda);  // initilizae the layerPopulation for the next layer && ACTIVE FLAGE
+                  layerDiscardPopulation[currentLayerIDX].CopyTo(layerPopulation[nextLayerIDX], 0);
+                  numberDiscard[currentLayerIDX] = 0;
                 }
                 else {
@@ -184 +234 @@
               }
               else {  // next layer is active
-                layerPopulation[i+1].CopyTo(layerJointPopulation[i + 1], 0);
-                layerJointPopulation[i].Where((x, idx) => idx == discardedIndividualIndex).ToArray().CopyTo(layerJointPopulation[i + 1], populationSize);
-                SmetricSelection(lambda, i + 1); // AGE and HVC 
-                numberDiscard[i] = 0;
+                layerPopulation[nextLayerIDX].CopyTo(layerJointPopulation[nextLayerIDX], 0);
+                layerJointPopulation[currentLayerIDX].Where((x, idx) => idx == discardedIndividualIndex).ToArray().CopyTo(layerJointPopulation[nextLayerIDX], populationSize);
+                SmetricSelection(lambda, nextLayerIDX); // AGE and HVC 
+                numberDiscard[currentLayerIDX] = 0;
               }
             }
-            layerPopulation[i].CopyTo(population, 0);  // BUG: should copy all the active layers to population.
+            layerPopulation[currentLayerIDX].CopyTo(population, 0);  // BUG: should copy all the active layers to population.
+              #endregion
+            }
           }
           else {
             // Some thing wrong? lol nothing wrong here ^_^
           }
-        }
-
+
+
+        }
+
+
+
+
+        int numberOfActiveLayer = activeLayer.Where(c => c).Count();
+        population = new IDynamicALPSSolution[populationSize * numberOfActiveLayer];
+        for (int i = 0; i < numberOfActiveLayer; i++) {
+          layerPopulation[i].CopyTo(population, i * populationSize);
+        }
 
 
@@ -203 +265 @@
         ExecuteOperation(executionContext, innerToken, analyze);
         // update Pareto-front approximation sets
-        UpdateParetoFronts();
+        // UpdateParetoFronts();
         // Show some results in the GUI.
         Results.AddOrUpdateResult("IdealPoint", new DoubleArray(IdealPoint));
@@ -209 +271 @@
         Results.AddOrUpdateResult("Evaluated Solutions", new IntValue(evaluatedSolutions));
 
-        var allLayerResults = new ResultCollection();
-        Results.AddOrUpdateResult("LayerResults", allLayerResults);
-        
+
+
+
+
+
+        // see if we already have a result collection for the layer results, and reuse that one
+        ResultCollection allLayerResults;
+        if (Results.TryGetValue("LayerResults", out IResult allLayerResultsResult)) {
+          allLayerResults = (ResultCollection)allLayerResultsResult.Value;
+        }
+        else {
+          allLayerResults = new ResultCollection();
+          Results.AddOrUpdateResult("LayerResults", allLayerResults);
+        }
+
+
+
+
+
         // run layer analyzers
         for (int i = 0; i < activeLayer.Length; ++i) {
@@ -217 +295 @@
           var scope = new Scope($"Layer {i}");
           var layer = layerPopulation[i];
-          var layerResults = new ResultCollection();
-          allLayerResults.Add(new Result(scope.Name, layerResults));
-
-          var problem = (MultiObjectiveTestFunctionProblem)Problem;
-          //scope.Variables.Add(new Variable("Individuals", new ItemArray<IScope>(layer.Select(x => (IScope)x.Individual))));
-          //scope.Variables.Add(new Variable("Qualities", new ItemArray<DoubleArray>(layer.Select(x => new DoubleArray(x.Qualities)))));
-          scope.SubScopes.AddRange(layer.Select(x => (IScope)x.Individual));
-          scope.Variables.Add(new Variable("Results", layerResults));
-          scope.Variables.Add(new Variable("TestFunction", problem.TestFunction));
-          var scatterPlotAnalyzer = new ScatterPlotAnalyzer();
-          scatterPlotAnalyzer.IndividualsParameter.ActualName = "RealVector";
-          var scattetPlotAnalyzerContext = executionContext.CreateChildOperation(scatterPlotAnalyzer, scope);
-          ExecuteOperation(executionContext, CancellationToken.None, scattetPlotAnalyzerContext);
-          scope.SubScopes.Clear();
-          scope.Variables.Clear();
+          var tmp = UpdateParetoFronts(layer, IdealPoint);
+
+          // update the results in a way that avoids creating a new result collection at each iteration
+          if (allLayerResults.TryGetValue(scope.Name, out IResult lRes)) {
+            var lr = (ResultCollection)lRes.Value;
+            foreach (var result in tmp) {
+              lr.AddOrUpdateResult(result.Name, result.Value);
+            }
+          }
+          else {
+            allLayerResults.AddOrUpdateResult(scope.Name, tmp);
+          }
+
+          var layerResults = (ResultCollection)allLayerResults[scope.Name].Value;
+
+          //var layerQualities = new ItemArray<DoubleArray>(layer.Select(x => new DoubleArray(x.Qualities)));
+          // var layerSolutions = new ItemArray<IItem>(layer.Select(x => (IItem)x.Individual.Clone()));
+
+          // only store the decision vectors
+          var layerSolutions = new ItemArray<IItem>(layer.Select(x => (IItem)((IScope)x.Individual).Variables["RealVector"].Value.Clone()));
+
+          var layerAges = new ItemArray<IntValue>(layer.Select(x => new IntValue(x.Age)));
+          //layerResults.AddOrUpdateResult("Objective values", layerQualities);
+          layerResults.AddOrUpdateResult("Decision vectors", layerSolutions);
+          layerResults.AddOrUpdateResult("Age", layerAges);
+
+          var tableObjectives = new DataTable("Objective values");
+          for (int j = 0; j < IdealPoint.Length; ++j) {
+            var row = new DataRow($"Objective {j + 1}");
+            row.Values.AddRange(layer.Select(x => x.Qualities[j]));
+            tableObjectives.Rows.Add(row);
+          }
+          layerResults.AddOrUpdateResult("Objective values", tableObjectives);
+
+
+          // historical HV 
+          DataTable hyperVolumeHistory;
+          if (layerResults.TryGetValue("Layer Hypervolume History", out IResult res)) {
+            hyperVolumeHistory = (DataTable)res.Value;
+          }
+          else {
+            hyperVolumeHistory = new DataTable("Layer Hypervolume History");
+            var hrow = new DataRow($"Layer {i}");
+            hrow.VisualProperties = new DataRowVisualProperties {
+              StartIndexZero = false,
+            };
+            hyperVolumeHistory.Rows.Add(hrow);
+            layerResults.AddOrUpdateResult("Layer Hypervolume History", hyperVolumeHistory);
+          }
+          //var front = layer.Select(x => x.Qualities);
+          var reference = ReferencePoint.ToArray();
+          //var hv = double.MinValue;
+
+          var layerQualities = layer.Select(x => x.Qualities).ToArray();
+          var layerPF = DominationCalculator<IDynamicALPSSolution>.CalculateBestParetoFront(layer, layerQualities, maximization);
+          var nondominatedLayer = NonDominatedSelect.GetDominatingVectors(layerPF.Select(x => x.Item2), reference, maximization, false);
+          var layerHV = nondominatedLayer.Any() ? Hypervolume.Calculate(nondominatedLayer, reference, maximization) : 0;
+          hyperVolumeHistory.Rows[$"Layer {i}"].Values.Add(layerHV);
+
+
+          // historical crossover probability 
+          DataTable crossoverProbabilityHistory;
+          if (layerResults.TryGetValue("CrossoverProbability History", out IResult resPm)) {
+            crossoverProbabilityHistory = (DataTable)resPm.Value;
+          }
+          else {
+            crossoverProbabilityHistory = new DataTable("CrossoverProbability History");
+            var hrowPm = new DataRow($"Layer {i}");
+            hrowPm.VisualProperties = new DataRowVisualProperties {
+              StartIndexZero = false,
+            };
+            crossoverProbabilityHistory.Rows.Add(hrowPm);
+            layerResults.AddOrUpdateResult("CrossoverProbability History", crossoverProbabilityHistory);
+          }
+          crossoverProbabilityHistory.Rows[$"Layer {i}"].Values.Add(layerCrossoverProbability[i]);
+
+
+
+          if (i == 1) {
+            DataTable wholeLayerHypervolumeHistrory;
+            var qualities = population.Select(x => x.Qualities).ToArray();
+            var pf = DominationCalculator<IDynamicALPSSolution>.CalculateBestParetoFront(population, qualities, maximization);
+            var nondominatedWhole = NonDominatedSelect.GetDominatingVectors(pf.Select(x => x.Item2), reference, maximization, false);
+            var hvWhole = nondominatedWhole.Any() ? Hypervolume.Calculate(nondominatedWhole, reference, maximization) : 0;
+            if (layerResults.TryGetValue("Hypervolume of the entire layers -- History", out IResult resHVWhole)) {
+              wholeLayerHypervolumeHistrory = (DataTable)resHVWhole.Value;
+            }
+            else {
+              wholeLayerHypervolumeHistrory = new DataTable("Hypervolume of the entire layers -- History");
+              var hrowWhole = new DataRow($"Layer {i}");
+              hrowWhole.VisualProperties = new DataRowVisualProperties {
+                StartIndexZero = false,
+              };
+              wholeLayerHypervolumeHistrory.Rows.Add(hrowWhole);
+              layerResults.AddOrUpdateResult("Hypervolume of the entire layers -- History", wholeLayerHypervolumeHistrory);
+            }
+            wholeLayerHypervolumeHistrory.Rows[$"Layer {i}"].Values.Add(hvWhole);
+          }
+          else {
+          }
+
         }
 
         // Update globalScope
         globalScope.SubScopes.Replace(population.Select(x => (IScope)x.Individual));
-
-
-        //// intilize the population for the next layer
-        //counterLayerALPS += 1;
-        //layerPopulation[counterLayerALPS] = new IDynamicALPSSolution[populationSize];
-        //layerPopulation[counterLayerALPS-1].CopyTo(layerPopulation[counterLayerALPS], 0); // DETELTE DUBGU 
-        //// BUG lambda should be different~~~~!!!!
-        //layerOffspringPopulation[counterLayerALPS] = new IDynamicALPSSolution[lambda];
       }
     }

branches/3057_DynamicALPS/HeuristicLab.Algorithms.DynamicALPS/3.4/DynamicALPSAlgorithmBase.cs

@@ -87 +87 @@
 
     [Storable]
+    protected double[] layerCrossoverProbability;
+
+    [Storable]
     protected IDynamicALPSSolution[][] layerDiscardPopulation;
+
+    [Storable]
+    protected IDynamicALPSSolution[] layerDiscardIndivdual;
+    
 
     [Storable]
@@ -116 +123 @@
     protected ExecutionState executionState;
 
-    private DoubleArray ReferencePoint {
+    protected DoubleArray ReferencePoint {
       get {
-        var problem = (MultiObjectiveTestFunctionProblem)Problem;
-        return problem.ReferencePoint;
+        if (Problem is MultiObjectiveTestFunctionProblem) {
+          var problem = (MultiObjectiveTestFunctionProblem)Problem;
+          return problem.ReferencePoint;
+        }
+        else {
+          return null;
+        }
       }
     }
@@ -136 +148 @@
     private const string RandomParameterName = "Random";
     private const string AnalyzerParameterName = "Analyzer";
+    
+    
     // MOEA-D parameters
     //private const string NeighbourSizeParameterName = "NeighbourSize";
@@ -147 +161 @@
     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";
 
 
@@ -199 +214 @@
     }
 
+    private IValueParameter<BoolValue> ALPSInitialzeLayerPopulationParameter {
+      get { return (IValueParameter<BoolValue>)Parameters[InitializeLayerPopulationMethodName]; }
+    }
+    
+
     private IValueParameter<IntValue> ResultPopulationSizeParameter {
       get { return (IValueParameter<IntValue>)Parameters[ResultPopulationSizeParameterName]; }
@@ -258 +278 @@
       get { return ALPSAgeGapParameter.Value; }
       set { ALPSAgeGapParameter.Value = value; }
+    }
+    public BoolValue UseAverageAge {
+      get { return ALPSInitialzeLayerPopulationParameter.Value; }
+      set { ALPSInitialzeLayerPopulationParameter.Value = value; }
     }
 
@@ -316 +340 @@
       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
@@ -532 +557 @@
       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);
@@ -550 +577 @@
         // 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 = ReferencePoint.ToArray();
+
+        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;
@@ -565 +605 @@
           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;
     }
 
@@ -581 +638 @@
       var maximumEvaluatedSolutions = MaximumEvaluatedSolutions.Value;
       var crossover = Crossover;
-      var crossoverProbability = CrossoverProbability.Value;
+      var crossoverProbability = layerCrossoverProbability[0];
       var mutator = Mutator;
       var mutationProbability = MutationProbability.Value;
@@ -604 +661 @@
       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
@@ -617 +687 @@
       else { // MUTATION   POLISHI
         if (childScope == null) {
-          offSpringAge = ages[mates[0]] + 1;
+          offSpringAge = ages[mates[0]];
         }
         else {
@@ -624 +694 @@
         var opMutation = executionContext.CreateChildOperation(mutator, childScope);
         ExecuteOperation(executionContext, innerToken, opMutation);
+        offSpringAge = offSpringAge + 1;
       }
       if (childScope != null) { // Evaluate the childScope
@@ -663 +734 @@
 
 
-
-    // Update the Pareto-front approximation set and scatter the solutions in PF approximation set. 
-    protected void UpdateParetoFronts() {
+      // 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 = population.Select(x => x.Qualities).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(population, qualities, maximization);
-
-      var pf2 = DominationCalculator<IDynamicALPSSolution>.CalculateAllParetoFronts(population, qualities, maximization, out int[] ranking);
+      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);
 
@@ -687 +759 @@
       //               True:  maximization problem 
       //               False: minimization problem
-      var reference = ReferencePoint.ToArray();
+      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;
@@ -696 +780 @@
 
       var elementNames = new List<string>() { "Pareto Front" };
-
-      ResultCollection results;
-      if (Results.ContainsKey("Hypervolume Analysis")) {
-        results = (ResultCollection)Results["Hypervolume Analysis"].Value;
+      var results = new ResultCollection();
+
+      ResultCollection innerResults;
+      if (results.ContainsKey("Hypervolume Analysis")) {
+        innerResults = (ResultCollection)results["Hypervolume Analysis"].Value;
       }
       else {
-        results = new ResultCollection();
-        Results.AddOrUpdateResult("Hypervolume Analysis", results);
+        innerResults = new ResultCollection();
+        results.AddOrUpdateResult("Hypervolume Analysis", innerResults);
       }
 
@@ -712 +797 @@
         if (error != OnlineCalculatorError.None) { r = double.NaN; }
         var resultName = "Pareto Front Analysis ";
-        if (!results.ContainsKey(resultName)) {
+        if (!innerResults.ContainsKey(resultName)) {
           sp = new ScatterPlot() {
-            VisualProperties = {
-              XAxisMinimumAuto = false, XAxisMinimumFixedValue = 0d, XAxisMaximumAuto = false, XAxisMaximumFixedValue = 1d,
-              YAxisMinimumAuto = false, YAxisMinimumFixedValue = 0d, YAxisMaximumAuto = false, YAxisMaximumFixedValue = 1d
-            }
+            //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 } });
-          results.AddOrUpdateResult(resultName, sp);
+          innerResults.AddOrUpdateResult(resultName, sp);
         }
         else {
-          sp = (ScatterPlot)results[resultName].Value;
+          sp = (ScatterPlot)innerResults[resultName].Value;
           sp.Rows[resultName].Points.Replace(points);
         }
@@ -751 +836 @@
           var pf2dPoints = pf2d.Select(x => new Point2D<double>(x.Item2[0], x.Item2[1]));
 
-          if (!results.ContainsKey(resultName)) {
+          if (!innerResults.ContainsKey(resultName)) {
             sp = new ScatterPlot() {
               VisualProperties = {
@@ -760 +845 @@
             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 } });
-            results.AddOrUpdateResult(resultName, sp);
+            innerResults.AddOrUpdateResult(resultName, sp);
           }
           else {
-            sp = (ScatterPlot)results[resultName].Value;
+            sp = (ScatterPlot)innerResults[resultName].Value;
             sp.Rows["Pareto Front"].Points.Replace(points);
             sp.Rows[$"Pareto Front [{c[0]}, {c[1]}]"].Points.Replace(pf2dPoints);
@@ -773 +858 @@
       }
       hypervolumes.RowNames = elementNames;
-      results.AddOrUpdateResult("Hypervolumes", hypervolumes);
+      innerResults.AddOrUpdateResult("Hypervolumes", hypervolumes);
+
+      return results;
     }