1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022018 Heuristic and Evolutionary Algorithms Laboratory (HEAL)


4  *


5  * This file is part of HeuristicLab.


6  *


7  * HeuristicLab is free software: you can redistribute it and/or modify


8  * it under the terms of the GNU General Public License as published by


9  * the Free Software Foundation, either version 3 of the License, or


10  * (at your option) any later version.


11  *


12  * HeuristicLab is distributed in the hope that it will be useful,


13  * but WITHOUT ANY WARRANTY; without even the implied warranty of


14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


15  * GNU General Public License for more details.


16  *


17  * You should have received a copy of the GNU General Public License


18  * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.


19  */


20  #endregion


21  using System;


22  using System.Collections.Generic;


23  using System.Linq;


24  using HeuristicLab.Common;


25 


26  namespace HeuristicLab.Optimization {


27  public static class HypervolumeCalculator {


28  public static double[] CalculateNadirPoint(IEnumerable<double[]> qualities, IReadOnlyList<bool> maximization) {


29  var res = maximization.Select(m => m ? double.MaxValue : double.MinValue).ToArray();


30  foreach (var quality in qualities)


31  for (var i = 0; i < quality.Length; i++)


32  if (maximization[i] == res[i] > quality[i])


33  res[i] = quality[i];


34  return res;


35  }


36 


37  /// <summary>


38  /// The Hypervolumemetric is defined as the HypervolumeCalculator enclosed between a given reference point,


39  /// that is fixed for every evaluation function and the evaluated qualities.


40  ///


41  /// Example:


42  /// r is the reference point at (11) and every point p is part of the evaluated qualities


43  /// The filled area labled HV is the 2 dimensional HypervolumeCalculator enclosed by this qualities.


44  ///


45  /// (01) (11)


46  /// + +r


47  ///  ###### HV ###


48  ///  p+######


49  ///  p+#####


50  ///  #####


51  ///  p+###


52  ///  p+


53  /// 


54  /// +1


55  /// (00) (10)


56  ///


57  /// Please note that in this example both dimensions are minimized. The reference point needs to be dominated by EVERY point in the evaluated qualities


58  ///


59  /// </summary>


60  ///


61  public static double CalculateHypervolume(IList<double[]> qualities, double[] referencePoint, IReadOnlyList<bool> maximization) {


62  qualities = qualities.Where(vec => DominationCalculator.Dominates(vec, referencePoint, maximization, false) == DominationResult.Dominates).ToArray();


63  if (qualities.Count== 0) return 0; //TODO computation for negative hypervolume?


64  if (maximization.Count == 2)


65  return Calculate2D(qualities, referencePoint, maximization);


66 


67  if (maximization.All(x => !x))


68  return CalculateMultiDimensional(qualities, referencePoint);


69  throw new NotImplementedException("HypervolumeCalculator calculation for more than two dimensions is supported only with minimization problems.");


70  }


71 


72 


73  /// <summary>


74  /// Caluclates the Hypervolume for a 2 dimensional problem


75  /// </summary>


76  /// <param name="front">All points within the front need to be NonDominated and need to dominate the reference point</param>


77  /// <param name="referencePoint"></param>


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


79  /// <returns></returns>


80  public static double Calculate2D(IList<double[]> front, double[] referencePoint, IReadOnlyList<bool> maximization) {


81  if (front == null) throw new ArgumentNullException("front");


82  if (referencePoint == null) throw new ArgumentNullException("referencePoint");


83  if (maximization == null) throw new ArgumentNullException("maximization");


84  if (!front.Any()) throw new ArgumentException("Front must not be empty.");


85  if (referencePoint.Length != 2) throw new ArgumentException("ReferencePoint must have exactly two dimensions.");


86 


87  var set = front.ToArray();


88  if (set.Any(s => s.Length != 2)) throw new ArgumentException("Points in qualities must have exactly two dimensions.");


89 


90  Array.Sort(set, new DimensionComparer(0, maximization[0]));


91 


92  double sum = 0;


93  for (var i = 0; i < set.Length  1; i++)


94  sum += Math.Abs(set[i][0]  set[i + 1][0]) * Math.Abs(set[i][1]  referencePoint[1]);


95  var lastPoint = set[set.Length  1];


96  sum += Math.Abs(lastPoint[0]  referencePoint[0]) * Math.Abs(lastPoint[1]  referencePoint[1]);


97 


98  return sum;


99  }


100 


101  public static double CalculateMultiDimensional(IList<double[]> front, double[] referencePoint) {


102  if (referencePoint == null  referencePoint.Length < 3) throw new ArgumentException("ReferencePoint unfit for complex HypervolumeCalculator calculation");


103 


104  var objectives = referencePoint.Length;


105  var fronList = front.ToList();


106  fronList.StableSort(new DimensionComparer(objectives  1, false));


107 


108  var regLow = Enumerable.Repeat(1E15, objectives).ToArray();


109  foreach (var p in fronList) {


110  for (var i = 0; i < regLow.Length; i++) {


111  if (p[i] < regLow[i]) regLow[i] = p[i];


112  }


113  }


114 


115  return Stream(regLow, referencePoint, fronList, 0, referencePoint[objectives  1], (int)Math.Sqrt(fronList.Count), objectives);


116  }


117 


118 


119  //within Stream a number of equality comparisons on double values are performed


120  //this is intentional and required


121  private static double Stream(double[] regionLow, double[] regionUp, List<double[]> front, int split, double cover, int sqrtNoPoints, int objectives) {


122  var coverOld = cover;


123  var coverIndex = 0;


124  var coverIndexOld = 1;


125  int c;


126  double result = 0;


127 


128  var dMeasure = GetMeasure(regionLow, regionUp, objectives);


129  while (cover == coverOld && coverIndex < front.Count()) {


130  if (coverIndexOld == coverIndex) break;


131  coverIndexOld = coverIndex;


132  if (Covers(front[coverIndex], regionLow, objectives)) {


133  cover = front[coverIndex][objectives  1];


134  result += dMeasure * (coverOld  cover);


135  } else coverIndex++;


136  }


137 


138  for (c = coverIndex; c > 0; c) if (front[c  1][objectives  1] == cover) coverIndex;


139  if (coverIndex == 0) return result;


140 


141  var allPiles = true;


142  var piles = new int[coverIndex];


143  for (var i = 0; i < coverIndex; i++) {


144  piles[i] = IsPile(front[i], regionLow, objectives);


145  if (piles[i] == 1) {


146  allPiles = false;


147  break;


148  }


149  }


150 


151  if (allPiles) {


152  var trellis = new double[regionUp.Length];


153  for (var j = 0; j < trellis.Length; j++) trellis[j] = regionUp[j];


154  double next;


155  var i = 0;


156  do {


157  var current = front[i][objectives  1];


158  do {


159  if (front[i][piles[i]] < trellis[piles[i]]) trellis[piles[i]] = front[i][piles[i]];


160  i++;


161  if (i < coverIndex) next = front[i][objectives  1];


162  else {


163  next = cover;


164  break;


165  }


166  }


167  while (next == current);


168  result += ComputeTrellis(regionLow, regionUp, trellis, objectives) * (next  current);


169  }


170  while (next != cover);


171  } else {


172  double bound = 1;


173  var boundaries = new double[coverIndex];


174  var noBoundaries = new double[coverIndex];


175  var boundIdx = 0;


176  var noBoundIdx = 0;


177 


178  do {


179  for (var i = 0; i < coverIndex; i++) {


180  var contained = ContainesBoundary(front[i], regionLow, split);


181  if (contained == 0) boundaries[boundIdx++] = front[i][split];


182  else if (contained == 1) noBoundaries[noBoundIdx++] = front[i][split];


183  }


184  if (boundIdx > 0) bound = GetMedian(boundaries, boundIdx);


185  else if (noBoundIdx > sqrtNoPoints) bound = GetMedian(noBoundaries, noBoundIdx);


186  else split++;


187  }


188  while (bound == 1.0);


189 


190  var pointsChildLow = new List<double[]>();


191  var pointsChildUp = new List<double[]>();


192  var regionUpC = new double[regionUp.Length];


193  for (var j = 0; j < regionUpC.Length; j++) regionUpC[j] = regionUp[j];


194  var regionLowC = new double[regionLow.Length];


195  for (var j = 0; j < regionLowC.Length; j++) regionLowC[j] = regionLow[j];


196 


197  for (var i = 0; i < coverIndex; i++) {


198  if (PartCovers(front[i], regionUpC, objectives)) pointsChildUp.Add(front[i]);


199  if (PartCovers(front[i], regionUp, objectives)) pointsChildLow.Add(front[i]);


200  }


201 


202  if (pointsChildUp.Count > 0) result += Stream(regionLow, regionUpC, pointsChildUp, split, cover, sqrtNoPoints, objectives);


203  if (pointsChildLow.Count > 0) result += Stream(regionLowC, regionUp, pointsChildLow, split, cover, sqrtNoPoints, objectives);


204  }


205  return result;


206  }


207 


208  private static double GetMedian(double[] vector, int length) {


209  return vector.Take(length).Median();


210  }


211 


212  private static double ComputeTrellis(double[] regionLow, double[] regionUp, double[] trellis, int objectives) {


213  var bs = new bool[objectives  1];


214  for (var i = 0; i < bs.Length; i++) bs[i] = true;


215 


216  double result = 0;


217  var noSummands = BinarayToInt(bs);


218  for (uint i = 1; i <= noSummands; i++) {


219  double summand = 1;


220  IntToBinary(i, bs);


221  var oneCounter = 0;


222  for (var j = 0; j < objectives  1; j++) {


223  if (bs[j]) {


224  summand *= regionUp[j]  trellis[j];


225  oneCounter++;


226  } else {


227  summand *= regionUp[j]  regionLow[j];


228  }


229  }


230  if (oneCounter % 2 == 0) result = summand;


231  else result += summand;


232  }


233  return result;


234  }


235 


236  private static void IntToBinary(uint i, bool[] bs) {


237  for (var j = 0; j < bs.Length; j++) bs[j] = false;


238  var rest = i;


239  var idx = 0;


240  while (rest != 0) {


241  bs[idx] = rest % 2 == 1;


242  rest = rest / 2;


243  idx++;


244  }


245  }


246 


247  private static uint BinarayToInt(bool[] bs) {


248  uint result = 0;


249  for (var i = 0; i < bs.Length; i++) {


250  result += bs[i] ? ((uint)1 << i) : 0;


251  }


252  return result;


253  }


254 


255  private static int IsPile(double[] cuboid, double[] regionLow, int objectives) {


256  var pile = cuboid.Length;


257  for (var i = 0; i < objectives  1; i++) {


258  if (cuboid[i] > regionLow[i]) {


259  if (pile != objectives) return 1;


260  pile = i;


261  }


262  }


263  return pile;


264  }


265 


266  private static double GetMeasure(double[] regionLow, double[] regionUp, int objectives) {


267  double volume = 1;


268  for (var i = 0; i < objectives  1; i++) {


269  volume *= (regionUp[i]  regionLow[i]);


270  }


271  return volume;


272  }


273 


274  private static int ContainesBoundary(double[] cub, double[] regionLow, int split) {


275  if (regionLow[split] >= cub[split]) return 1;


276  else {


277  for (var j = 0; j < split; j++) {


278  if (regionLow[j] < cub[j]) return 1;


279  }


280  }


281  return 0;


282  }


283 


284  private static bool PartCovers(double[] v, double[] regionUp, int objectives) {


285  for (var i = 0; i < objectives  1; i++) {


286  if (v[i] >= regionUp[i]) return false;


287  }


288  return true;


289  }


290 


291  private static bool Covers(double[] v, double[] regionLow, int objectives) {


292  for (var i = 0; i < objectives  1; i++) {


293  if (v[i] > regionLow[i]) return false;


294  }


295  return true;


296  }


297 


298  private class DimensionComparer : IComparer<double[]> {


299  private readonly int dimension;


300  private readonly int descending;


301 


302  public DimensionComparer(int dimension, bool descending) {


303  this.dimension = dimension;


304  this.descending = descending ? 1 : 1;


305  }


306 


307  public int Compare(double[] x, double[] y) {


308  return x[dimension].CompareTo(y[dimension]) * descending;


309  }


310  }


311  }


312  } 
