1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022016 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.Problems.MultiObjectiveTestFunctions {


27 


28  public static class Hypervolume {


29 


30  /// <summary>


31  /// The Hyprevolumemetric is defined as the Hypervolume enclosed between a given reference point,


32  /// that is fixed for every evaluation function and the evaluated front.


33  ///


34  /// Example:


35  /// r is the reference Point at (11) and every Point p is part of the evaluated front


36  /// The filled Area labled HV is the 2 dimensional Hypervolume enclosed by this front.


37  ///


38  /// (01) (11)


39  /// + +r


40  ///  ###### HV ###


41  ///  p+######


42  ///  p+#####


43  ///  #####


44  ///  p+###


45  ///  p+


46  /// 


47  /// +1


48  /// (00) (10)


49  ///


50  /// Please note that in this example both dimensions are minimized. The reference Point need to be dominated by EVERY point in the evaluated front


51  ///


52  /// </summary>


53  ///


54  public static double Calculate(IEnumerable<double[]> front, double[] referencePoint, bool[] maximization) {


55  front = NonDominatedSelect.GetDominatingVectors(front, referencePoint, maximization, false);


56  if (!front.Any()) throw new ArgumentException("No point in the front dominates the referencePoint");


57  if (maximization.Length == 2)


58  return Calculate2D(front, referencePoint, maximization);


59 


60  if (Array.TrueForAll(maximization, x => !x))


61  return CalculateMulitDimensional(front, referencePoint);


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


63  }


64 


65 


66  private static double Calculate2D(IEnumerable<double[]> front, double[] referencePoint, bool[] maximization) {


67  if (front == null) throw new ArgumentNullException("Front must not be null.");


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


69 


70  if (referencePoint == null) throw new ArgumentNullException("ReferencePoint must not be null.");


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


72 


73  double[][] set = front.ToArray();


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


75 


76  Array.Sort<double[]>(set, new Utilities.DimensionComparer(0, maximization[0]));


77 


78  double sum = 0;


79  for (int i = 0; i < set.Length  1; i++) {


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


81  }


82 


83  double[] lastPoint = set[set.Length  1];


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


85 


86  return sum;


87  }


88 


89 


90 


91  private static double CalculateMulitDimensional(IEnumerable<double[]> front, double[] referencePoint) {


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


93 


94  int objectives = referencePoint.Length;


95  var fronList = front.ToList();


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


97 


98  double[] regLow = Enumerable.Repeat(1E15, objectives).ToArray();


99  foreach (double[] p in fronList) {


100  for (int i = 0; i < regLow.Length; i++) {


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


102  }


103  }


104 


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


106  }


107 


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


109  double coverOld = cover;


110  int coverIndex = 0;


111  int coverIndexOld = 1;


112  int c;


113  double result = 0;


114 


115  double dMeasure = GetMeasure(regionLow, regionUp, objectives);


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


117  if (coverIndexOld == coverIndex) break;


118  coverIndexOld = coverIndex;


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


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


121  result += dMeasure * (coverOld  cover);


122  } else coverIndex++;


123 


124  }


125 


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


127  if (coverIndex == 0) return result;


128 


129  bool allPiles = true;


130  int[] piles = new int[coverIndex];


131  for (int i = 0; i < coverIndex; i++) {


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


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


134  allPiles = false;


135  break;


136  }


137  }


138 


139  if (allPiles) {


140  double[] trellis = new double[regionUp.Length];


141  for (int j = 0; j < trellis.Length; j++) trellis[j] = regionUp[j];


142  double current = 0;


143  double next = 0;


144  int i = 0;


145  do {


146  current = front[i][objectives  1];


147  do {


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


149  i++;


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


151  else { next = cover; break; }


152  } while (next == current);


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


154  } while (next != cover);


155  } else {


156  double bound = 1;


157  double[] boundaries = new double[coverIndex];


158  double[] noBoundaries = new double[coverIndex];


159  int boundIdx = 0;


160  int noBoundIdx = 0;


161 


162  do {


163  for (int i = 0; i < coverIndex; i++) {


164  int contained = ContainesBoundary(front[i], regionLow, split);


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


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


167  }


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


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


170  else split++;


171  } while (bound == 1.0);


172 


173  List<double[]> pointsChildLow, pointsChildUp;


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


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


176  double[] regionUpC = new double[regionUp.Length];


177  for (int j = 0; j < regionUpC.Length; j++) regionUpC[j] = regionUp[j];


178  double[] regionLowC = new double[regionLow.Length];


179  for (int j = 0; j < regionLowC.Length; j++) regionLowC[j] = regionLow[j];


180 


181  for (int i = 0; i < coverIndex; i++) {


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


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


184  }


185  //this could/should be done in Parallel


186 


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


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


189  }


190  return result;


191  }


192 


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


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


195  }


196 


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


198  bool[] bs = new bool[objectives  1];


199  for (int i = 0; i < bs.Length; i++) bs[i] = true;


200 


201  double result = 0;


202  uint noSummands = BinarayToInt(bs);


203  int oneCounter; double summand;


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


205  summand = 1;


206  IntToBinary(i, bs);


207  oneCounter = 0;


208  for (int j = 0; j < objectives  1; j++) {


209  if (bs[j]) {


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


211  oneCounter++;


212  } else {


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


214  }


215  }


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


217  else result += summand;


218 


219  }


220  return result;


221  }


222 


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


224  for (int j = 0; j < bs.Length; j++) bs[j] = false;


225  uint rest = i;


226  int idx = 0;


227  while (rest != 0) {


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


229  rest = rest / 2;


230  idx++;


231  }


232 


233  }


234 


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


236  uint result = 0;


237  for (int i = 0; i < bs.Length; i++) {


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


239  }


240  return result;


241  }


242 


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


244  int pile = cuboid.Length;


245  for (int i = 0; i < objectives  1; i++) {


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


247  if (pile != objectives) return 1;


248  pile = i;


249  }


250  }


251  return pile;


252  }


253 


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


255  double volume = 1;


256  for (int i = 0; i < objectives  1; i++) {


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


258  }


259  return volume;


260  }


261 


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


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


264  else {


265  for (int j = 0; j < split; j++) {


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


267  }


268  }


269  return 0;


270  }


271 


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


273  for (int i = 0; i < objectives  1; i++) {


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


275  }


276  return true;


277  }


278 


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


280  for (int i = 0; i < objectives  1; i++) {


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


282  }


283  return true;


284  }


285 


286 


287  }


288  }


289 

