Free cookie consent management tool by TermsFeed Policy Generator

source: branches/2521_ProblemRefactoring/HeuristicLab.Optimization/3.3/MultiObjective/HypervolumeCalculator.cs @ 17309

Last change on this file since 17309 was 17309, checked in by abeham, 4 years ago

#2521: Refactored maximization property for multi-objective problems

File size: 12.1 KB
Line 
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
21using System;
22using System.Collections.Generic;
23using System.Linq;
24using HeuristicLab.Common;
25
26namespace 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 Hypervolume-metric 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 (1|1) 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    /// (0|1)                (1|1)
46    ///   +      +-------------r
47    ///   |      |###### HV ###|
48    ///   |      p------+######|
49    ///   |             p+#####|
50    ///   |              |#####|
51    ///   |              p-+###|
52    ///   |                p---+
53    ///   |                 
54    ///   +--------------------1
55    /// (0|0)                (1|0)
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 Non-Dominated 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}
Note: See TracBrowser for help on using the repository browser.