[16171] | 1 | #region License Information
|
---|
| 2 | /* HeuristicLab
|
---|
| 3 | * Copyright (C) 2002-2018 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 |
|
---|
| 28 | public static class HypervolumeCalculator {
|
---|
| 29 |
|
---|
| 30 | public static double[] CalculateNadirPoint(IEnumerable<double[]> qualities, bool[] maximization){
|
---|
| 31 | var res = maximization.Select(m => m ? double.MaxValue : double.MinValue).ToArray();
|
---|
| 32 | foreach (var quality in qualities)
|
---|
| 33 | for (var i = 0; i < quality.Length; i++)
|
---|
| 34 | if (maximization[i] == res[i] > quality[i])
|
---|
| 35 | res[i] = quality[i];
|
---|
| 36 | return res;
|
---|
| 37 | }
|
---|
| 38 |
|
---|
| 39 | /// <summary>
|
---|
| 40 | /// The Hyprevolume-metric is defined as the HypervolumeCalculator enclosed between a given reference point,
|
---|
| 41 | /// that is fixed for every evaluation function and the evaluated qualities.
|
---|
| 42 | ///
|
---|
| 43 | /// Example:
|
---|
| 44 | /// r is the reference point at (1|1) and every point p is part of the evaluated qualities
|
---|
| 45 | /// The filled area labled HV is the 2 dimensional HypervolumeCalculator enclosed by this qualities.
|
---|
| 46 | ///
|
---|
| 47 | /// (0|1) (1|1)
|
---|
| 48 | /// + +-------------r
|
---|
| 49 | /// | |###### HV ###|
|
---|
| 50 | /// | p------+######|
|
---|
| 51 | /// | p+#####|
|
---|
| 52 | /// | |#####|
|
---|
| 53 | /// | p-+###|
|
---|
| 54 | /// | p---+
|
---|
| 55 | /// |
|
---|
| 56 | /// +--------------------1
|
---|
| 57 | /// (0|0) (1|0)
|
---|
| 58 | ///
|
---|
| 59 | /// Please note that in this example both dimensions are minimized. The reference point needs to be dominated by EVERY point in the evaluated qualities
|
---|
| 60 | ///
|
---|
| 61 | /// </summary>
|
---|
| 62 | ///
|
---|
| 63 | public static double CalculateHypervolume(double[][] qualities, double[] referencePoint, bool[] maximization){
|
---|
| 64 | qualities = qualities.Where(vec => DominationCalculator.Dominates(vec, referencePoint, maximization, false) == DominationResult.Dominates).ToArray();
|
---|
[16310] | 65 | if (qualities.Length == 0 ) return 0;//TODO computation for negative hypervolume?
|
---|
[16171] | 66 | if (maximization.Length == 2)
|
---|
| 67 | return Calculate2D(qualities, referencePoint, maximization);
|
---|
| 68 |
|
---|
| 69 | if (Array.TrueForAll(maximization, x => !x))
|
---|
| 70 | return CalculateMultiDimensional(qualities, referencePoint);
|
---|
| 71 | throw new NotImplementedException("HypervolumeCalculator calculation for more than two dimensions is supported only with minimization problems.");
|
---|
| 72 | }
|
---|
| 73 |
|
---|
[16310] | 74 |
|
---|
| 75 | /// <summary>
|
---|
| 76 | /// Caluclates the Hypervolume for a 2 dimensional problem
|
---|
| 77 | /// </summary>
|
---|
| 78 | /// <param name="front">All points within the front need to be Non-Dominated and need to dominate the reference point</param>
|
---|
| 79 | /// <param name="referencePoint"></param>
|
---|
| 80 | /// <param name="maximization"></param>
|
---|
| 81 | /// <returns></returns>
|
---|
[16171] | 82 | public static double Calculate2D(double[][] front, double[] referencePoint, bool[] maximization) {
|
---|
[16310] | 83 | if (front == null) throw new ArgumentNullException(nameof(front));
|
---|
| 84 | if (referencePoint == null) throw new ArgumentNullException(nameof(referencePoint));
|
---|
| 85 | if (maximization == null) throw new ArgumentNullException(nameof(maximization));
|
---|
[16171] | 86 | if (!front.Any()) throw new ArgumentException("Front must not be empty.");
|
---|
| 87 | if (referencePoint.Length != 2) throw new ArgumentException("ReferencePoint must have exactly two dimensions.");
|
---|
| 88 |
|
---|
| 89 | var set = front.ToArray();
|
---|
| 90 | if (set.Any(s => s.Length != 2)) throw new ArgumentException("Points in qualities must have exactly two dimensions.");
|
---|
| 91 |
|
---|
| 92 | Array.Sort(set, new DimensionComparer(0, maximization[0]));
|
---|
| 93 |
|
---|
| 94 | double sum = 0;
|
---|
| 95 | for (var i = 0; i < set.Length - 1; i++)
|
---|
[16310] | 96 | sum += Math.Abs(set[i][0] - set[i + 1][0]) * Math.Abs(set[i][1] - referencePoint[1]);
|
---|
[16171] | 97 | var lastPoint = set[set.Length - 1];
|
---|
| 98 | sum += Math.Abs(lastPoint[0] - referencePoint[0]) * Math.Abs(lastPoint[1] - referencePoint[1]);
|
---|
| 99 |
|
---|
| 100 | return sum;
|
---|
| 101 | }
|
---|
| 102 |
|
---|
| 103 | public static double CalculateMultiDimensional(double[][] front, double[] referencePoint) {
|
---|
| 104 | if (referencePoint == null || referencePoint.Length < 3) throw new ArgumentException("ReferencePoint unfit for complex HypervolumeCalculator calculation");
|
---|
| 105 |
|
---|
| 106 | var objectives = referencePoint.Length;
|
---|
| 107 | var fronList = front.ToList();
|
---|
| 108 | fronList.StableSort(new DimensionComparer(objectives - 1, false));
|
---|
| 109 |
|
---|
| 110 | var regLow = Enumerable.Repeat(1E15, objectives).ToArray();
|
---|
| 111 | foreach (var p in fronList) {
|
---|
| 112 | for (var i = 0; i < regLow.Length; i++) {
|
---|
| 113 | if (p[i] < regLow[i]) regLow[i] = p[i];
|
---|
| 114 | }
|
---|
| 115 | }
|
---|
| 116 |
|
---|
| 117 | return Stream(regLow, referencePoint, fronList, 0, referencePoint[objectives - 1], (int)Math.Sqrt(fronList.Count), objectives);
|
---|
| 118 | }
|
---|
| 119 |
|
---|
[16310] | 120 |
|
---|
| 121 | //within Stream a number of equality comparisons on double values are performed
|
---|
| 122 | //this is intentional and required
|
---|
[16171] | 123 | private static double Stream(double[] regionLow, double[] regionUp, List<double[]> front, int split, double cover, int sqrtNoPoints, int objectives) {
|
---|
| 124 | var coverOld = cover;
|
---|
| 125 | var coverIndex = 0;
|
---|
| 126 | var coverIndexOld = -1;
|
---|
| 127 | int c;
|
---|
| 128 | double result = 0;
|
---|
| 129 |
|
---|
| 130 | var dMeasure = GetMeasure(regionLow, regionUp, objectives);
|
---|
| 131 | while (cover == coverOld && coverIndex < front.Count()) {
|
---|
| 132 | if (coverIndexOld == coverIndex) break;
|
---|
| 133 | coverIndexOld = coverIndex;
|
---|
| 134 | if (Covers(front[coverIndex], regionLow, objectives)) {
|
---|
| 135 | cover = front[coverIndex][objectives - 1];
|
---|
| 136 | result += dMeasure * (coverOld - cover);
|
---|
| 137 | } else coverIndex++;
|
---|
| 138 | }
|
---|
| 139 |
|
---|
| 140 | for (c = coverIndex; c > 0; c--) if (front[c - 1][objectives - 1] == cover) coverIndex--;
|
---|
| 141 | if (coverIndex == 0) return result;
|
---|
| 142 |
|
---|
| 143 | var allPiles = true;
|
---|
| 144 | var piles = new int[coverIndex];
|
---|
| 145 | for (var i = 0; i < coverIndex; i++) {
|
---|
[16310] | 146 | piles[i] = IsPile(front[i], regionLow, objectives);
|
---|
[16171] | 147 | if (piles[i] == -1) {
|
---|
| 148 | allPiles = false;
|
---|
| 149 | break;
|
---|
| 150 | }
|
---|
| 151 | }
|
---|
| 152 |
|
---|
| 153 | if (allPiles) {
|
---|
| 154 | var trellis = new double[regionUp.Length];
|
---|
| 155 | for (var j = 0; j < trellis.Length; j++) trellis[j] = regionUp[j];
|
---|
[16310] | 156 | double next;
|
---|
[16171] | 157 | var i = 0;
|
---|
| 158 | do {
|
---|
[16310] | 159 | var current = front[i][objectives - 1];
|
---|
[16171] | 160 | do {
|
---|
| 161 | if (front[i][piles[i]] < trellis[piles[i]]) trellis[piles[i]] = front[i][piles[i]];
|
---|
| 162 | i++;
|
---|
| 163 | if (i < coverIndex) next = front[i][objectives - 1];
|
---|
| 164 | else { next = cover; break; }
|
---|
| 165 | } while (next == current);
|
---|
| 166 | result += ComputeTrellis(regionLow, regionUp, trellis, objectives) * (next - current);
|
---|
| 167 | } while (next != cover);
|
---|
| 168 | } else {
|
---|
| 169 | double bound = -1;
|
---|
| 170 | var boundaries = new double[coverIndex];
|
---|
| 171 | var noBoundaries = new double[coverIndex];
|
---|
| 172 | var boundIdx = 0;
|
---|
| 173 | var noBoundIdx = 0;
|
---|
| 174 |
|
---|
| 175 | do {
|
---|
| 176 | for (var i = 0; i < coverIndex; i++) {
|
---|
| 177 | var contained = ContainesBoundary(front[i], regionLow, split);
|
---|
| 178 | if (contained == 0) boundaries[boundIdx++] = front[i][split];
|
---|
| 179 | else if (contained == 1) noBoundaries[noBoundIdx++] = front[i][split];
|
---|
| 180 | }
|
---|
| 181 | if (boundIdx > 0) bound = GetMedian(boundaries, boundIdx);
|
---|
| 182 | else if (noBoundIdx > sqrtNoPoints) bound = GetMedian(noBoundaries, noBoundIdx);
|
---|
| 183 | else split++;
|
---|
| 184 | } while (bound == -1.0);
|
---|
| 185 |
|
---|
[16310] | 186 | var pointsChildLow = new List<double[]>();
|
---|
| 187 | var pointsChildUp = new List<double[]>();
|
---|
[16171] | 188 | var regionUpC = new double[regionUp.Length];
|
---|
| 189 | for (var j = 0; j < regionUpC.Length; j++) regionUpC[j] = regionUp[j];
|
---|
| 190 | var regionLowC = new double[regionLow.Length];
|
---|
| 191 | for (var j = 0; j < regionLowC.Length; j++) regionLowC[j] = regionLow[j];
|
---|
| 192 |
|
---|
| 193 | for (var i = 0; i < coverIndex; i++) {
|
---|
| 194 | if (PartCovers(front[i], regionUpC, objectives)) pointsChildUp.Add(front[i]);
|
---|
| 195 | if (PartCovers(front[i], regionUp, objectives)) pointsChildLow.Add(front[i]);
|
---|
| 196 | }
|
---|
| 197 |
|
---|
| 198 | if (pointsChildUp.Count > 0) result += Stream(regionLow, regionUpC, pointsChildUp, split, cover, sqrtNoPoints, objectives);
|
---|
| 199 | if (pointsChildLow.Count > 0) result += Stream(regionLowC, regionUp, pointsChildLow, split, cover, sqrtNoPoints, objectives);
|
---|
| 200 | }
|
---|
| 201 | return result;
|
---|
| 202 | }
|
---|
| 203 |
|
---|
| 204 | private static double GetMedian(double[] vector, int length) {
|
---|
| 205 | return vector.Take(length).Median();
|
---|
| 206 | }
|
---|
| 207 |
|
---|
| 208 | private static double ComputeTrellis(double[] regionLow, double[] regionUp, double[] trellis, int objectives) {
|
---|
| 209 | var bs = new bool[objectives - 1];
|
---|
| 210 | for (var i = 0; i < bs.Length; i++) bs[i] = true;
|
---|
| 211 |
|
---|
| 212 | double result = 0;
|
---|
| 213 | var noSummands = BinarayToInt(bs);
|
---|
| 214 | for (uint i = 1; i <= noSummands; i++) {
|
---|
[16310] | 215 | double summand = 1;
|
---|
[16171] | 216 | IntToBinary(i, bs);
|
---|
[16310] | 217 | var oneCounter = 0;
|
---|
[16171] | 218 | for (var j = 0; j < objectives - 1; j++) {
|
---|
| 219 | if (bs[j]) {
|
---|
| 220 | summand *= regionUp[j] - trellis[j];
|
---|
| 221 | oneCounter++;
|
---|
| 222 | } else {
|
---|
| 223 | summand *= regionUp[j] - regionLow[j];
|
---|
| 224 | }
|
---|
| 225 | }
|
---|
| 226 | if (oneCounter % 2 == 0) result -= summand;
|
---|
| 227 | else result += summand;
|
---|
| 228 |
|
---|
| 229 | }
|
---|
| 230 | return result;
|
---|
| 231 | }
|
---|
| 232 |
|
---|
| 233 | private static void IntToBinary(uint i, bool[] bs) {
|
---|
| 234 | for (var j = 0; j < bs.Length; j++) bs[j] = false;
|
---|
| 235 | var rest = i;
|
---|
| 236 | var idx = 0;
|
---|
| 237 | while (rest != 0) {
|
---|
| 238 | bs[idx] = rest % 2 == 1;
|
---|
| 239 | rest = rest / 2;
|
---|
| 240 | idx++;
|
---|
| 241 | }
|
---|
| 242 | }
|
---|
| 243 |
|
---|
| 244 | private static uint BinarayToInt(bool[] bs) {
|
---|
| 245 | uint result = 0;
|
---|
| 246 | for (var i = 0; i < bs.Length; i++) {
|
---|
| 247 | result += bs[i] ? ((uint)1 << i) : 0;
|
---|
| 248 | }
|
---|
| 249 | return result;
|
---|
| 250 | }
|
---|
| 251 |
|
---|
[16310] | 252 | private static int IsPile(double[] cuboid, double[] regionLow, int objectives) {
|
---|
[16171] | 253 | var pile = cuboid.Length;
|
---|
| 254 | for (var i = 0; i < objectives - 1; i++) {
|
---|
| 255 | if (cuboid[i] > regionLow[i]) {
|
---|
| 256 | if (pile != objectives) return 1;
|
---|
| 257 | pile = i;
|
---|
| 258 | }
|
---|
| 259 | }
|
---|
| 260 | return pile;
|
---|
| 261 | }
|
---|
| 262 |
|
---|
| 263 | private static double GetMeasure(double[] regionLow, double[] regionUp, int objectives) {
|
---|
| 264 | double volume = 1;
|
---|
| 265 | for (var i = 0; i < objectives - 1; i++) {
|
---|
| 266 | volume *= (regionUp[i] - regionLow[i]);
|
---|
| 267 | }
|
---|
| 268 | return volume;
|
---|
| 269 | }
|
---|
| 270 |
|
---|
| 271 | private static int ContainesBoundary(double[] cub, double[] regionLow, int split) {
|
---|
| 272 | if (regionLow[split] >= cub[split]) return -1;
|
---|
| 273 | else {
|
---|
| 274 | for (var j = 0; j < split; j++) {
|
---|
| 275 | if (regionLow[j] < cub[j]) return 1;
|
---|
| 276 | }
|
---|
| 277 | }
|
---|
| 278 | return 0;
|
---|
| 279 | }
|
---|
| 280 |
|
---|
| 281 | private static bool PartCovers(double[] v, double[] regionUp, int objectives) {
|
---|
| 282 | for (var i = 0; i < objectives - 1; i++) {
|
---|
| 283 | if (v[i] >= regionUp[i]) return false;
|
---|
| 284 | }
|
---|
| 285 | return true;
|
---|
| 286 | }
|
---|
| 287 |
|
---|
| 288 | private static bool Covers(double[] v, double[] regionLow, int objectives) {
|
---|
| 289 | for (var i = 0; i < objectives - 1; i++) {
|
---|
| 290 | if (v[i] > regionLow[i]) return false;
|
---|
| 291 | }
|
---|
| 292 | return true;
|
---|
| 293 | }
|
---|
| 294 |
|
---|
| 295 | private class DimensionComparer : IComparer<double[]> {
|
---|
| 296 | private readonly int dimension;
|
---|
| 297 | private readonly int descending;
|
---|
| 298 |
|
---|
| 299 | public DimensionComparer(int dimension, bool descending) {
|
---|
| 300 | this.dimension = dimension;
|
---|
| 301 | this.descending = descending ? -1 : 1;
|
---|
| 302 | }
|
---|
| 303 |
|
---|
| 304 | public int Compare(double[] x, double[] y){
|
---|
| 305 | return x[dimension].CompareTo(y[dimension]) * descending;
|
---|
| 306 | }
|
---|
| 307 | }
|
---|
| 308 |
|
---|
| 309 | }
|
---|
| 310 | }
|
---|
| 311 |
|
---|