1  using System;


2  using System.Collections.Generic;


3  using System.Linq;


4  using HeuristicLab.Core;


5  using HeuristicLab.Optimization;


6 


7  namespace HeuristicLab.Algorithms.NSGA3


8  {


9  public static class NSGA3Selection


10  {


11  private const double EPSILON = 10e6; // a tiny number that is greater than 0


12 


13  /// <summary>


14  /// TODO: Description If the fitnesses of all solutions in <paramref name="rt" /> do not


15  /// have the same number of objectives, undefined behavior happens.


16  /// </summary>


17  /// <param name="rt"></param>


18  /// <param name="referencePoints"></param>


19  /// <param name="maximization"></param>


20  /// <returns></returns>


21  public static List<Solution> SelectSolutionsForNextGeneration(List<Solution> rt, List<ReferencePoint> referencePoints, bool[] maximization, IRandom random)


22  {


23  int N = rt.Count / 2; // number of solutions in a generation


24  int numberOfObjectives = rt.First().Fitness.Length;


25  List<Solution> st = new List<Solution>();


26 


27  // Do nondominated sort


28  var qualities = Utility.ToFitnessMatrix(rt);


29  // compute the pareto fronts using the DominationCalculator and discard the qualities


30  // part in the inner tuples


31  var fronts = DominationCalculator<Solution>.CalculateAllParetoFronts(rt.ToArray(), qualities, maximization, out int[] rank, true)


32  .Select(list => new List<Solution>(list.Select(pair => pair.Item1))).ToList();


33 


34  int lf = 1; // last front to be included


35  for (int i = 0; i < fronts.Count && st.Count < N; i++)


36  {


37  st = Utility.Concat(st, fronts[i]);


38  lf = i;


39  }


40  // remove useless fronts


41  fronts.RemoveRange(lf + 1, fronts.Count  lf  1);


42 


43  List<Solution> nextGeneration = new List<Solution>();


44  if (st.Count == N) // no special selection needs to be done


45  nextGeneration = st;


46  else


47  {


48  // Add every front to the next generation list except for the one that is added only partially


49  nextGeneration = new List<Solution>();


50  for (int f = 0; f < lf; f++)


51  nextGeneration = Utility.Concat(nextGeneration, fronts[f]);


52  int k = N  nextGeneration.Count;


53  Normalize(st, numberOfObjectives);


54  Associate(fronts, referencePoints);


55  List<Solution> solutionsToAdd = Niching(k, referencePoints, random);


56  nextGeneration = Utility.Concat(nextGeneration, solutionsToAdd);


57  }


58 


59  return nextGeneration;


60  }


61 


62  private static void Normalize(List<Solution> st, int numberOfObjectives)


63  {


64  // Find the ideal point


65  double[] idealPoint = new double[numberOfObjectives];


66  for (int j = 0; j < numberOfObjectives; j++)


67  {


68  // Compute ideal point


69  idealPoint[j] = Utility.Min(s => s.Fitness[j], st);


70 


71  // Translate objectives


72  foreach (var solution in st)


73  solution.Fitness[j] = idealPoint[j];


74  }


75 


76  // Find the extreme points


77  Solution[] extremePoints = new Solution[numberOfObjectives];


78  for (int j = 0; j < numberOfObjectives; j++)


79  {


80  // Compute extreme points


81  double[] weights = new double[numberOfObjectives];


82  for (int i = 0; i < numberOfObjectives; i++) weights[i] = EPSILON;


83  weights[j] = 1;


84 


85  extremePoints[j] = Utility.ArgMin(s => ASF(s.Fitness, weights), st);


86  }


87 


88  // Compute intercepts


89  List<double> intercepts = GetIntercepts(extremePoints.ToList());


90 


91  // Normalize objectives


92  NormalizeObjectives(st, intercepts, idealPoint);


93  }


94 


95  private static void NormalizeObjectives(List<Solution> st, List<double> intercepts, double[] idealPoint)


96  {


97  int numberOfObjectives = st.First().Fitness.Length;


98 


99  foreach (var solution in st)


100  {


101  for (int i = 0; i < numberOfObjectives; i++)


102  {


103  if (Math.Abs(intercepts[i]  idealPoint[i]) > EPSILON)


104  {


105  solution.Fitness[i] = solution.Fitness[i] / (intercepts[i]  idealPoint[i]);


106  }


107  else


108  {


109  solution.Fitness[i] = solution.Fitness[i] / EPSILON;


110  }


111  }


112  }


113  }


114 


115  private static void Associate(List<List<Solution>> fronts, List<ReferencePoint> referencePoints)


116  {


117  for (int f = 0; f < fronts.Count; f++)


118  {


119  foreach (var solution in fronts[f])


120  {


121  // find reference point for which the perpendicular distance to the current


122  // solution is the lowest


123  var rpAndDist = Utility.MinArgMin(rp => GetPerpendicularDistance(rp.Values, solution.Fitness), referencePoints);


124  // associated reference point


125  var arp = rpAndDist.Item1;


126  // distance to that reference point


127  var dist = rpAndDist.Item2;


128 


129  if (f + 1 != fronts.Count)


130  {


131  // Todo: Add member for reference point on index min_rp


132  arp.NumberOfAssociatedSolutions++;


133  }


134  else


135  {


136  // Todo: Add potential member for reference point on index min_rp


137  arp.AddPotentialAssociatedSolution(solution, dist);


138  }


139  }


140  }


141  }


142 


143  private static List<Solution> Niching(int k, List<ReferencePoint> referencePoints, IRandom random)


144  {


145  List<Solution> solutions = new List<Solution>();


146  while (solutions.Count != k)


147  {


148  ReferencePoint min_rp = FindNicheReferencePoint(referencePoints, random);


149 


150  Solution chosen = SelectClusterMember(min_rp);


151  if (chosen == null)


152  {


153  referencePoints.Remove(min_rp);


154  }


155  else


156  {


157  min_rp.NumberOfAssociatedSolutions++;


158  min_rp.RemovePotentialAssociatedSolution(chosen);


159  solutions.Add(chosen);


160  }


161  }


162 


163  return solutions;


164  }


165 


166  private static ReferencePoint FindNicheReferencePoint(List<ReferencePoint> referencePoints, IRandom random)


167  {


168  // the minimum number of associated solutions for a reference point over all reference points


169  int minNumber = Utility.Min(rp => rp.NumberOfAssociatedSolutions, referencePoints);


170 


171  // the reference points that share the number of associated solutions where that number


172  // is equal to minNumber


173  List<ReferencePoint> minAssociatedReferencePoints = new List<ReferencePoint>();


174  foreach (var referencePoint in referencePoints)


175  if (referencePoint.NumberOfAssociatedSolutions == minNumber)


176  minAssociatedReferencePoints.Add(referencePoint);


177 


178  if (minAssociatedReferencePoints.Count > 1)


179  return minAssociatedReferencePoints[random.Next(minAssociatedReferencePoints.Count)];


180  else


181  return minAssociatedReferencePoints.Single();


182  }


183 


184  private static Solution SelectClusterMember(ReferencePoint referencePoint)


185  {


186  Solution chosen = null;


187  if (referencePoint.HasPotentialMember())


188  {


189  if (referencePoint.NumberOfAssociatedSolutions == 0)


190  chosen = referencePoint.FindClosestMember();


191  else


192  chosen = referencePoint.RandomMember();


193  }


194  return chosen;


195  }


196 


197  private static double GetPerpendicularDistance(double[] values, double[] fitness)


198  {


199  double numerator = 0;


200  double denominator = 0;


201  for (int i = 0; i < values.Length; i++)


202  {


203  numerator += values[i] * fitness[i];


204  denominator += Math.Pow(values[i], 2);


205  }


206  double k = numerator / denominator;


207 


208  double d = 0;


209  for (int i = 0; i < values.Length; i++)


210  {


211  d += Math.Pow(k * values[i]  fitness[i], 2);


212  }


213  return Math.Sqrt(d);


214  }


215 


216  private static double ASF(double[] x, double[] weight)


217  {


218  int numberOfObjectives = x.Length;


219  List<int> dimensions = new List<int>();


220  for (int i = 0; i < numberOfObjectives; i++)


221  dimensions.Add(i);


222 


223  return Utility.Max((dim) => x[dim] / weight[dim], dimensions);


224  }


225 


226  private static List<double> GetIntercepts(List<Solution> extremePoints)


227  {


228  int numberOfObjectives = extremePoints.First().Fitness.Length;


229 


230  // Check whether there are duplicate extreme points. This might happen but the original


231  // paper does not mention how to deal with it.


232  bool duplicate = false;


233  for (int i = 0; !duplicate && i < extremePoints.Count; i++)


234  {


235  for (int j = i + 1; !duplicate && j < extremePoints.Count; j++)


236  {


237  // maybe todo: override Equals method of solution?


238  duplicate = extremePoints[i].Equals(extremePoints[j]);


239  }


240  }


241 


242  List<double> intercepts = new List<double>();


243 


244  if (duplicate)


245  { // cannot construct the unique hyperplane (this is a casual method to deal with the condition)


246  for (int f = 0; f < numberOfObjectives; f++)


247  {


248  // extreme_points[f] stands for the individual with the largest value of


249  // objective f


250  intercepts.Add(extremePoints[f].Fitness[f]);


251  }


252  }


253  else


254  {


255  // Find the equation of the hyperplane


256  List<double> b = new List<double>(); //(pop[0].objs().size(), 1.0);


257  for (int i = 0; i < numberOfObjectives; i++)


258  {


259  b.Add(1.0);


260  }


261 


262  List<List<double>> a = new List<List<double>>();


263  foreach (Solution s in extremePoints)


264  {


265  List<double> aux = new List<double>();


266  for (int i = 0; i < numberOfObjectives; i++)


267  aux.Add(s.Fitness[i]);


268  a.Add(aux);


269  }


270  List<double> x = GaussianElimination(a, b);


271 


272  // Find intercepts


273  for (int f = 0; f < numberOfObjectives; f++)


274  {


275  intercepts.Add(1.0 / x[f]);


276  }


277  }


278 


279  return intercepts;


280  }


281 


282  private static List<double> GaussianElimination(List<List<double>> a, List<double> b)


283  {


284  List<double> x = new List<double>();


285 


286  int n = a.Count;


287  for (int i = 0; i < n; i++)


288  a[i].Add(b[i]);


289 


290  for (int @base = 0; @base < n  1; @base++)


291  for (int target = @base + 1; target < n; target++)


292  {


293  double ratio = a[target][@base] / a[@base][@base];


294  for (int term = 0; term < a[@base].Count; term++)


295  a[target][term] = a[target][term]  a[@base][term] * ratio;


296  }


297 


298  for (int i = 0; i < n; i++)


299  x.Add(0.0);


300 


301  for (int i = n  1; i >= 0; i)


302  {


303  for (int known = i + 1; known < n; known++)


304  a[i][n] = a[i][n]  a[i][known] * x[known];


305  x[i] = a[i][n] / a[i][i];


306  }


307 


308  return x;


309  }


310  }


311  } 
