source: branches/HeuristicLab.Problems.MultiObjectiveTestFunctions/HeuristicLab.Problems.MultiObjectiveTestFunctions/3.3/Calculators/MultiDimensionalHypervolume.cs @ 14018

#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 <http://www.gnu.org/licenses/>.
 */
#endregion

using System;
using System.Collections.Generic;
using System.Linq;
using HeuristicLab.Common;

namespace HeuristicLab.Problems.MultiObjectiveTestFunctions {

  public class MultiDimensionalHypervolume {
    /// <summary>
    /// 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 
    /// 
    /// This is an efficient calculation for a multidimensional Hypervolume as described in 
    /// "Faster S-Metric Calculation by Considering Dominated
    /// Hypervolume as Klee’s Measure Problem" by
    /// Nicola Beume and Günter Rudolph.
    /// 
    /// 
    /// A reference Impelementation in C++ can be found at http://image.diku.dk/shark/doxygen_pages/html/_hypervolume_calculator_8h_source.html
    /// </summary>
    public static double Calculate(IEnumerable<double[]> points, IEnumerable<double> reference) {
      double[] referencePoint = reference.ToArray<double>();
      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<double[]> lpoints = new List<double[]>();
      foreach (double[] p in points) {
        lpoints.Add(p);
      }
      lpoints.StableSort(new Utilities.DimensionComparer(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<double[]> 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<double[]> 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<double[]> pointsChildLow, pointsChildUp;
        pointsChildLow = new List<double[]>();
        pointsChildUp = new List<double[]>();
        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;
    }


  }
}