source: branches/2943_MOBasicProblem_MOCMAES/HeuristicLab.Optimization/3.3/MultiObjective/HypervolumeCalculator.cs @ 16171

Last change on this file since 16171 was 16171, checked in by bwerth, 2 years ago

#2943 worked on MOBasicProblem - added Interfaces;reworked MOCalculators; several minor changes

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