source: branches/2520_PersistenceReintegration/HeuristicLab.Problems.TestFunctions.MultiObjective/3.3/Calculators/HyperVolume.cs @ 16453

Last change on this file since 16453 was 16453, checked in by jkarder, 2 years ago

#2520: updated year of copyrights

File size: 10.9 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2019 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.Problems.TestFunctions.MultiObjective {
27
28  public static class Hypervolume {
29
30    /// <summary>
31    /// The Hyprevolume-metric 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 (1|1) 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    /// (0|1)                (1|1)
39    ///   +      +-------------r
40    ///   |      |###### HV ###|
41    ///   |      p------+######|
42    ///   |             p+#####|
43    ///   |              |#####|
44    ///   |              p-+###|
45    ///   |                p---+
46    ///   |                 
47    ///   +--------------------1
48    /// (0|0)                (1|0)
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
Note: See TracBrowser for help on using the repository browser.