source: branches/HeuristicLab.Problems.MultiObjectiveTestFunctions/HeuristicLab.Problems.MultiObjectiveTestFunctions/3.3/Comparators/HyperVolumeFast.cs @ 13620

using System;
using System.Collections.Generic;
using HeuristicLab.Common;
using HeuristicLab.Problems.MultiObjectiveTestFunctions.Comparators;

namespace HeuristicLab.Problems.MultiObjectiveTestFunctions {
  /// <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 
  /// 
  /// </summary>
  public class FastHypervolume {

    private double[] reference;
    private bool[] maximization;
    public FastHypervolume(double[] reference, bool[] maximization) {
      if (reference.Length != 2) throw new NotSupportedException("Only 2-dimensional cases are supported yet");
      this.reference = reference;
      this.maximization = maximization;
    }

    public double Compare(IEnumerable<double[]> front, IEnumerable<double[]> optimalFront, double[,] bounds) {
      return GetHypervolume(optimalFront, bounds) - GetHypervolume(front, bounds);

    }


    /// <summary>
    /// 
    /// </summary>
    /// <param name="front"></param>
    /// <param name="bounds">also called region</param>
    /// <returns></returns>
    public double GetHypervolume(IEnumerable<double[]> front, double[,] bounds) {
      //TODO what to do if set contains dominated points
      if (front == null) throw new ArgumentException("Fronts must not be null");
      var list = new List<double[]>();
      list.AddRange(front);
      list.Sort(Utilities.getDimensionComparer(reference.Length - 1, maximization[reference.Length - 1]));
      var results = new ExPrivates(this);
      GetHV(SpanRegion(), list, 0, reference[reference.Length - 1], results);
      return results.Volume;

    }

    private double[][] SpanRegion(double[,] bounds, int length) {
      if (bounds.GetLength(1) != 2) throw new ArgumentException();
      double[][] copy = new double[length][];
      for (int i = 0; i < copy.Length; i++) {
        copy[i] = new double[2];
        for (int j = 0; j < 2; j++) {
          copy[i][j] = bounds[i%bounds.GetLength(0),j];
        }
      }
      return copy;
    }

    private double[][] SpanRegion() {
      double[][] copy = new double[reference.Length][];
      for (int i = 0; i < copy.Length; i++) {
        copy[i] = new double[2];
        copy[i][0] = 0;
        copy[i][1] = reference[i];
      }
      return copy;

    }


    private double[][] DeepClone(double[][] array) {
      double[][] copy = new double[array.Length][];
      for (int i = 0; i < copy.Length; i++) {
        copy[i] = new double[array[i].Length];
        for (int j = 0; j < array[i].Length; j++) {
          copy[i][j] = array[i][j];
        }
      }
      return copy;
    }

    private class ExPrivates {
      private double volume;
      public double Volume {
        get { return volume; }
        set { volume = value; }
      }
      private int d;

      public int D {
        get { return d; }
      }
      public ExPrivates(FastHypervolume hv) {
        volume = 0;
        d = hv.reference.Length;
      }



    }

    private void GetHV(double[][] region, List<double[]> points, int split, double cover, ExPrivates results) {
      double coverNew = cover;
      int coverIndex = 0;
      bool allPiles = true;
      double bound = -1;

      /* is the region completely covered? */
      while (coverNew == cover && coverIndex != points.Count) {
        if (covers(points[coverIndex], region, results)) {
          coverNew = points[coverIndex][results.D];
          results.Volume += GetMeasure(region,results) * (cover - coverNew);
        } else coverIndex++;
      }
      if (coverIndex == 0) return;

      for (int i = 0; i < coverIndex; i++) { if (checkPile(points[i], region, results) == -1) allPiles = false; }

      if (allPiles) {
        /* calculate volume by sweeping along dimension d */
        var trellis = new double[reference.Length];
        int i = 0;
        for (int j = 0; j < results.D-1; j++) trellis[j] = reference[j];

        double next;//bernhard
        do {
          double current = points[i][results.D-1];
          do {
            int pile = checkPile(points[i], region,results);
            if (points[i][pile] < trellis[pile]) trellis[pile] = points[i][pile];
            i++;
            if (i < coverIndex - 1) next = points[i][results.D-1]; else next = coverNew;
          } while (current == next);
          results.Volume += measure(trellis, region, results) * (next - current);
        } while (next != coverNew);

      } else {
        /* split region in two children regions */
        do {
          var intersect = new List<Double>(); var nonIntersect = new List<Double>();
          for (int i = 0; i < coverIndex; i++) {
            var intersection = intersects(points[i], region, split);
            if (intersection == 1) intersect.Add(points[i][split]);
            if (intersection == 0) nonIntersect.Add(points[i][split]);
          }
          if (intersect.Count != 0) bound = intersect.Median();
          else if (nonIntersect.Count > Math.Sqrt(points.Count)) bound = nonIntersect.Median();
          else split++;
        } while (bound == -1);
      }


      var regionC = DeepClone(region);
      regionC[split][1] = bound;
      var pointsC = new List<double[]>();
      for (int i = 0; i < coverIndex; i++) {
        if (partCovers(points[i], regionC, results)) move(points,i, pointsC);
      }
      if (pointsC.Count != 0) GetHV(regionC, pointsC, split, coverNew, results);
      reinsert(pointsC, points);
      regionC = region;
      regionC[split][0] = bound;
      pointsC = new List<double[]>();
      for (int i = 1; i < coverIndex; i++) {
        if (partCovers(points[i], regionC, results)) move(points,i, pointsC);
      }
      if (pointsC.Count != 0) GetHV(regionC, pointsC, split, coverNew, results);
      reinsert(pointsC, points);
    }

    private void reinsert(List<double[]> pointsC, List<double[]> points) {
      points.AddRange(pointsC);
      pointsC.Clear();
    }

    private double GetMeasure(double[][] region,ExPrivates result) {
      double volume = 1.0;
      // for ( std::size_t i = 0; i < regionLow.size(); i++ ) {
      for (int i = 1; i < result.D; i++) {
      volume *= (region[1][i] - region[0][i]);
      }
      return (volume);
    }

    private void move(List<double[]> points, int i, List<double[]> pointsC) {
      double[] v = points[i];
      points.Remove(v);
      pointsC.Add(v);
    }

    private double measure(double[] trellis, double[][] region, ExPrivates results) {
      double volume=0;
      bool[] indicator = new bool[results.D];
      for (int i = 0; i < results.D-1; i++) indicator[i] = true;
      int numberSummands = integerValue(indicator);

      for (int i = 0; i <= numberSummands; i++) {
        indicator = binaryValue(i);
        int oneCounter = 0;
        double summand = 0;
        for (int j = 1; j < results.D; j++) {
          if (indicator[i] == true) {
            summand += region[1][j] - trellis[j];
            oneCounter++;
          } else {
            summand += region[1][j] - region[0][j];
          }
        }
        if(oneCounter%2 == 0) {
          volume -= summand;
        } else {
          volume += summand;
        }
      }
      return volume;
    }

    private int integerValue(bool[] binary) {
      int sum=0;
      foreach (bool b in binary) {
        sum = (sum << 1) + (b ? 1 : 0);
      }
      return sum;
    }

    private bool[] binaryValue(int integer) {
      bool[] res = new bool[32];
      int i = 0;
      while (integer != 0) {
        res[i++]= (integer & 1)==1;
        integer >>= 1;
      }
      return res;

    }

    private int checkPile(double[] point, double[][] region, ExPrivates results) {
      int pile = -1;
      for (int j = 0; j < results.D-1; j++) {
        if (point[j] > region[j][0]) {
          if (pile != -1) return -1;
        }
        pile = j;
      }
      return pile;

    }

    private int intersects(double[] point, double[][] region, int split) {
      if (region[split][0] >= point[split]) return -1;
      for (int j = 1; j < split; j++) {
        if (point[j] > region[j][0]) return 1;
      }
      return 0;

    }

    private bool covers(double[] point, double[][] region, ExPrivates results) {
      for (int j = 1; j < results.D; j++) {
        if (point[j] > region[j][0]) return false;
      }
      return true;

    }

    private bool partCovers(double[] point, double[][] region, ExPrivates results) {
      for (int j = 1; j < results.D; j++) {
        if (point[j] >= region[j][1]) return false;
      }
      return true;
    }

    public static double GetHypervolume(IEnumerable<double[]> front, double[] reference, bool[] maximization, double[,] bounds) {
      FastHypervolume comp = new FastHypervolume(reference, maximization);
      return comp.GetHypervolume(front,bounds);
    }

    public static double GetDistance(IEnumerable<double[]> front, IEnumerable<double[]> optimalFront, double[] reference, bool[] maximization, double[,] bounds) {
      return GetHypervolume(optimalFront, reference, maximization, bounds) - GetHypervolume(front, reference, maximization,bounds);
    }

    private void CheckConsistency(double[] point, int dim) {
      if (!maximization[dim] && point[dim] > reference[dim]) throw new Exception("Reference Point must be dominated by all points of the front");
      if (maximization[dim] && point[dim] < reference[dim]) throw new Exception("Reference Point must be dominated by all points of the front");
    }







   


  }
}