#region License Information
/* HeuristicLab
* Copyright (C) 2002-2016 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 .
*/
#endregion
using HeuristicLab.Common;
using System;
using System.Collections.Generic;
using System.Linq;
namespace HeuristicLab.Problems.MultiObjectiveTestFunctions
{
public class Hypervolume
{
///
/// The Hyprevolume-metric is defined as the Hypervolume enclosed between a given reference point,
/// that is fixed for every evaluation function and the evaluated front.
///
/// Example:
/// r is the reference Point at (1|1) and every Point p is part of the evaluated front
/// The filled Area labled HV is the 2 diensional Hypervolume enclosed by this front.
///
/// (0|1) (1|1)
/// + +-------------r
/// | |###### HV ###|
/// | p------+######|
/// | p+#####|
/// | |#####|
/// | p-+###|
/// | p---+
/// |
/// +--------------------1
/// (0|0) (1|0)
///
/// Please note that in this example both dimensions are minimized. The reference Point need to be dominated by EVERY point in the evaluated front
///
///
///
public static double Calculate(IEnumerable points, IEnumerable reference, bool[] maximization)
{
var front = NonDominatedSelect.removeNonReferenceDominatingVectors(points, reference.ToArray(), maximization, false);
if (maximization.Length == 2)
{ //Hypervolume analysis only with 2 objectives for now
return Hypervolume.Calculate2D(front, reference, maximization);
}
else if (Array.TrueForAll(maximization, x => !x))
{
return Hypervolume.CalculateMD(front, reference);
}
else
{
throw new NotImplementedException("Hypervolume calculation for more than two dimensions is supported only with minimization problems");
}
}
private static double Calculate2D(IEnumerable front, IEnumerable reference, bool[] maximization)
{
List list = new List();
foreach (double d in reference) list.Add(d);
double[] refp = list.ToArray();
if (front == null) throw new ArgumentException("Fronts must not be null");
double[][] set = front.ToArray(); //Still no Good
if (set.Length == 0) throw new ArgumentException("Fronts must not be empty");
if (refp.Length != set[0].Length) throw new ArgumentException("Front and referencepoint need to be of the same dimensionality");
Array.Sort(set, Utilities.getDimensionComparer(0, maximization[0]));
double[] last = set[set.Length - 1];
CheckConsistency(last, 0, refp, maximization);
CheckConsistency(last, 1, refp, maximization);
double sum = 0;
for (int i = 0; i < set.Length - 1; i++)
{
CheckConsistency(set[i], 1, refp, maximization);
sum += Math.Abs((set[i][0] - set[i + 1][0])) * Math.Abs((set[i][1] - refp[1]));
}
sum += Math.Abs(refp[0] - last[0]) * Math.Abs(refp[1] - last[1]);
return sum;
}
private static void CheckConsistency(double[] point, int dim, double[] reference, bool[] maximization)
{
if (!maximization[dim] && point[dim] > reference[dim]) throw new ArgumentException("Reference Point must be dominated by all points of the front");
if (maximization[dim] && point[dim] < reference[dim]) throw new ArgumentException("Reference Point must be dominated by all points of the front");
if (point.Length != 2) throw new ArgumentException("Only 2-dimensional cases are supported yet");
}
private static double CalculateMD(IEnumerable points, IEnumerable reference)
{
double[] referencePoint = reference.ToArray();
if (referencePoint == null || referencePoint.Length < 3) throw new ArgumentException("ReferencePoint unfit for complex Hypervolume calculation");
if (!IsDominated(referencePoint, points))
{
throw new ArgumentException("ReferencePoint unfit for complex Hypervolume calculation");
}
int objectives = referencePoint.Length;
List lpoints = new List();
foreach (double[] p in points)
{
lpoints.Add(p);
}
lpoints.StableSort(Utilities.getDimensionComparer(objectives - 1, false));
double[] regLow = new double[objectives];
for (int i = 0; i < objectives; i++)
{
regLow[i] = 1E15;
}
foreach (double[] p in lpoints)
{
for (int i = 0; i < regLow.Length; i++)
{
if (regLow[i] > p[i]) regLow[i] = p[i];
}
}
return Stream(regLow, referencePoint, lpoints, 0, referencePoint[objectives - 1], (int)Math.Sqrt(points.Count()), objectives);
}
private static bool IsDominated(double[] referencePoint, IEnumerable points)
{
foreach (double[] point in points)
{
for (int i = 0; i < referencePoint.Length; i++)
{
if (referencePoint[i] < point[i])
{
return false;
}
}
}
return true;
}
private static double Stream(double[] regionLow, double[] regionUp, List points, int split, double cover, int sqrtNoPoints, int objectives)
{
double coverOld = cover;
int coverIndex = 0;
int coverIndexOld = -1;
int c;
double result = 0;
double dMeasure = GetMeasure(regionLow, regionUp, objectives);
while (cover == coverOld && coverIndex < points.Count())
{
if (coverIndexOld == coverIndex) break;
coverIndexOld = coverIndex;
if (Covers(points[coverIndex], regionLow, objectives))
{
cover = points[coverIndex][objectives - 1];
result += dMeasure * (coverOld - cover);
}
else coverIndex++;
}
for (c = coverIndex; c > 0; c--) if (points[c - 1][objectives - 1] == cover) coverIndex--;
if (coverIndex == 0) return result;
bool allPiles = true;
int[] piles = new int[coverIndex];
for (int i = 0; i < coverIndex; i++)
{
piles[i] = IsPile(points[i], regionLow, regionUp, objectives);
if (piles[i] == -1)
{
allPiles = false;
break;
}
}
if (allPiles)
{
double[] trellis = new double[regionUp.Length];
for (int j = 0; j < trellis.Length; j++) trellis[j] = regionUp[j];
double current = 0;
double next = 0;
int i = 0;
do
{
current = points[i][objectives - 1];
do
{
if (points[i][piles[i]] < trellis[piles[i]]) trellis[piles[i]] = points[i][piles[i]];
i++;
if (i < coverIndex) next = points[i][objectives - 1];
else { next = cover; break; }
} while (next == current);
result += ComputeTrellis(regionLow, regionUp, trellis, objectives) * (next - current);
} while (next != cover);
}
else
{
double bound = -1;
double[] boundaries = new double[coverIndex];
double[] noBoundaries = new double[coverIndex];
int boundIdx = 0;
int noBoundIdx = 0;
do
{
for (int i = 0; i < coverIndex; i++)
{
int contained = ContainesBoundary(points[i], regionLow, split);
if (contained == 0) boundaries[boundIdx++] = points[i][split];
else if (contained == 1) noBoundaries[noBoundIdx++] = points[i][split];
}
if (boundIdx > 0) bound = GetMedian(boundaries, boundIdx);
else if (noBoundIdx > sqrtNoPoints) bound = GetMedian(noBoundaries, noBoundIdx);
else split++;
} while (bound == -1.0);
List pointsChildLow, pointsChildUp;
pointsChildLow = new List();
pointsChildUp = new List();
double[] regionUpC = new double[regionUp.Length];
for (int j = 0; j < regionUpC.Length; j++) regionUpC[j] = regionUp[j];
double[] regionLowC = new double[regionLow.Length];
for (int j = 0; j < regionLowC.Length; j++) regionLowC[j] = regionLow[j];
for (int i = 0; i < coverIndex; i++)
{
if (PartCovers(points[i], regionUpC, objectives)) pointsChildUp.Add(points[i]);
if (PartCovers(points[i], regionUp, objectives)) pointsChildLow.Add(points[i]);
}
//this could/should be done in Parallel
if (pointsChildUp.Count() > 0) result += Stream(regionLow, regionUpC, pointsChildUp, split, cover, sqrtNoPoints, objectives);
if (pointsChildLow.Count() > 0) result += Stream(regionLowC, regionUp, pointsChildLow, split, cover, sqrtNoPoints, objectives);
}
return result;
}
private static double GetMedian(double[] vector, int length)
{
if (vector.Length != length)
{
double[] vec = new double[length];
Array.Copy(vector, vec, length);
vector = vec;
}
return vector.Median();
}
private static double ComputeTrellis(double[] regionLow, double[] regionUp, double[] trellis, int objectives)
{
bool[] bs = new bool[objectives - 1];
for (int i = 0; i < bs.Length; i++) bs[i] = true;
double result = 0;
uint noSummands = BinarayToInt(bs);
int oneCounter; double summand;
for (uint i = 1; i <= noSummands; i++)
{
summand = 1;
IntToBinary(i, bs);
oneCounter = 0;
for (int j = 0; j < objectives - 1; j++)
{
if (bs[j])
{
summand *= regionUp[j] - trellis[j];
oneCounter++;
}
else
{
summand *= regionUp[j] - regionLow[j];
}
}
if (oneCounter % 2 == 0) result -= summand;
else result += summand;
}
return result;
}
private static void IntToBinary(uint i, bool[] bs)
{
for (int j = 0; j < bs.Length; j++) bs[j] = false;
uint rest = i;
int idx = 0;
while (rest != 0)
{
bs[idx] = rest % 2 == 1;
rest = rest / 2;
idx++;
}
}
private static uint BinarayToInt(bool[] bs)
{
uint result = 0;
for (int i = 0; i < bs.Length; i++)
{
result += bs[i] ? ((uint)1 << i) : 0;
}
return result;
}
private static int IsPile(double[] cuboid, double[] regionLow, double[] regionUp, int objectives)
{
int pile = cuboid.Length;
for (int i = 0; i < objectives - 1; i++)
{
if (cuboid[i] > regionLow[i])
{
if (pile != objectives) return 1;
pile = i;
}
}
return pile;
}
private static double GetMeasure(double[] regionLow, double[] regionUp, int objectives)
{
double volume = 1;
for (int i = 0; i < objectives - 1; i++)
{
volume *= (regionUp[i] - regionLow[i]);
}
return volume;
}
private static int ContainesBoundary(double[] cub, double[] regionLow, int split)
{
if (regionLow[split] >= cub[split]) return -1;
else
{
for (int j = 0; j < split; j++)
{
if (regionLow[j] < cub[j]) return 1;
}
}
return 0;
}
private static bool PartCovers(double[] v, double[] regionUp, int objectives)
{
for (int i = 0; i < objectives - 1; i++)
{
if (v[i] >= regionUp[i]) return false;
}
return true;
}
private static bool Covers(double[] v, double[] regionLow, int objectives)
{
for (int i = 0; i < objectives - 1; i++)
{
if (v[i] > regionLow[i]) return false;
}
return true;
}
}
}