source: branches/HeuristicLab.Analysis.AlgorithmBehavior/HeuristicLab.Analysis.AlgorithmBehavior.Analyzers/3.3/ConvexHullMeasures.cs @ 10287

#region License Information
/* HeuristicLab
 * Copyright (C) 2002-2013 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 Mehroz;
using MIConvexHull;

namespace HeuristicLab.Analysis.AlgorithmBehavior.Analyzers {
  public static class ConvexHullMeasures {
    //Calculates the volumne of a d-simplex
    public static double CalculateSimplexVolume(List<double[]> simplex) {
      int dim = simplex.First().Length;
      double[,] diffs = new double[dim, dim];

      for (int i = 1; i < simplex.Count; i++) {
        for (int j = 0; j < dim; j++) {
          diffs[j, i - 1] = simplex[i][j] - simplex[0][j];
        }
      }

      double det = Math.Abs(alglib.rmatrixdet(diffs));
      double result = det / dim.Fact();

      return result;
    }

    public static Fraction CalculateSimplexVolumeBig(List<double[]> simplex) {
      int dim = simplex.First().Length;
      double[,] diffs = new double[dim, dim];

      for (int i = 1; i < simplex.Count; i++) {
        for (int j = 0; j < dim; j++) {
          diffs[j, i - 1] = simplex[i][j] - simplex[0][j];
        }
      }

      //calculate determinante from lu decomposition
      int[] pivots = null;
      alglib.rmatrixlu(ref diffs, dim, dim, out pivots);

      Fraction det = 1;
      for (int i = 0; i < dim; i++) {
        det = det * diffs[i, i];
      }

      if (det < 0)
        det = -1 * det;

      Fraction result = det / dim.Fact();
      return result;
    }

    //calculates the volume of a convex d-polytope
    //decomposition based on delaunay triangulation
    public static double CalculateVolume(List<double[]> convexHull) {
      double volume = 0.0;

      int dim = convexHull.First().Length;
      if (dim > convexHull.Count)
        throw new ArgumentException("Nr. of points for volume calculation must be greater than dimension", "convexHull");

      ITriangulation<DefaultVertex, DefaultTriangulationCell<DefaultVertex>> triangulation = null;
      try {
        //Under certain circumstances MIConvexHull is not able to calculate the triangulation
        triangulation = Triangulation.CreateDelaunay(ConvertToVertex(convexHull));
      }
      catch {
        return double.NaN;
      }
      var simplices = ConvertFromTriangulation(triangulation);
      foreach (var simplex in simplices) {
        volume += CalculateSimplexVolume(simplex.ToList());
      }

      return volume;
    }

    public static Fraction CalculateVolumeBig(List<double[]> convexHull) {
      Fraction volume = new Fraction(0.0);

      int dim = convexHull.First().Length;
      if (dim > convexHull.Count)
        throw new ArgumentException("Nr. of points for volume calculation must be greater than dimension", "convexHull");

      ITriangulation<DefaultVertex, DefaultTriangulationCell<DefaultVertex>> triangulation = null;
      try {
        //Under certain circumstances MIConvexHull is not able to calculate the triangulation
        triangulation = Triangulation.CreateDelaunay(ConvertToVertex(convexHull));
      }
      catch {
        return Fraction.NaN;
      }
      var simplices = ConvertFromTriangulation(triangulation);
      foreach (var simplex in simplices) {
        volume += CalculateSimplexVolumeBig(simplex.ToList());
      }

      return volume;
    }

    private static List<double[][]> ConvertFromTriangulation(ITriangulation<DefaultVertex, DefaultTriangulationCell<DefaultVertex>> triangulation) {
      List<double[][]> results = new List<double[][]>();

      foreach (var cell in triangulation.Cells) {
        results.Add(cell.Vertices.Select(x => x.Position).ToArray());
      }

      return results;
    }

    private static List<DefaultVertex> ConvertToVertex(List<double[]> data) {
      List<DefaultVertex> result = new List<DefaultVertex>();
      for (int i = 0; i < data.Count(); i++) {
        double[] d = data[i];
        DefaultVertex vertex = new DefaultVertex();
        vertex.Position = d;
        result.Add(vertex);
      }
      return result;
    }

    public static double[] CalculateCentroids(List<double[]> convexHull) {
      int n = convexHull.Count;
      int dim = convexHull.First().Length;
      double[] centroid = new double[dim];

      for (int i = 0; i < dim; i++) {
        for (int j = 0; j < n; j++) {
          centroid[i] += convexHull[j][i];
        }
        centroid[i] /= n;
      }
      return centroid;
    }

    //measures to calculate: AVG, SUM, STD.DEV
    public static double[] CalculateMovementDistances(List<double[]> centroids) {
      double[] distances = new double[centroids.Count - 1];

      // it is assumed that the centroids are sorted chronological
      for (int i = 0; i < centroids.Count - 1; i++) {
        distances[i] = centroids[i].EuclideanDistance(centroids[i + 1]);
      }

      return distances;
    }

    public static double CalculateOverallMovementDistances(List<double[]> centroids) {
      // it is assumed that the centroids are sorted chronological
      return centroids[0].EuclideanDistance(centroids[centroids.Count - 1]);
    }

    //calculate theta from two vectors (dot product)
    public static double CalculateTheta(double[] a, double[] b) {
      return Math.Acos(a.DotProduct(b) / (a.EuclideanNorm() * b.EuclideanNorm()));
    }

    //calculate theta using law of cosines
    public static double CalculateTheta(double[] a, double[] b, double[] c) {
      double ab = a.EuclideanDistance(b);
      double bc = b.EuclideanDistance(c);
      double ca = c.EuclideanDistance(a);

      //return degrees
      return Math.Acos((ab * ab + bc * bc - ca * ca) / (2 * ab * bc)) * 180.0 / Math.PI;
    }

    //measures to calculate: AVG, SUM, STD.DEV
    public static double[] CalculateCentroidsMotion(List<double[]> centroids) {
      double[] thetas = new double[centroids.Count - 2];

      for (int i = 1; i < centroids.Count - 1; i++) {
        thetas[i - 1] = CalculateTheta(centroids[i - 1], centroids[i], centroids[i + 1]);
      }

      return thetas;
    }

    public static double CalculateMaxDiameter(List<double[]> convexHull) {
      double maxDist = double.MinValue;

      for (int i = 0; i < convexHull.Count; i++) {
        for (int j = 0; j < convexHull.Count; j++) {
          if (i != j) {
            double dist = convexHull[i].EuclideanDistance(convexHull[j]);
            if (dist > maxDist)
              maxDist = dist;
          }
        }
      }
      return maxDist;
    }

    private static List<double[]> GenerateHyperCube(double[] bounds, int dim) {
      if (bounds.Length != 2)
        throw new ArgumentException("Only bounds of length 2 are supported!", "bounds");

      int numPoints = (int)Math.Pow(2, dim);

      var hyperCube = new List<double[]>(numPoints);
      var binaryVectors = new List<string>(numPoints);
      for (int i = 0; i < numPoints; i++) {
        string binary = Convert.ToString(i, 2);
        binaryVectors.Add(binary);
      }

      for (int i = 0; i < numPoints; i++) {
        for (int j = 0; j < dim - binaryVectors[i].Length; j++) {
          //TODO: do this properly
          binaryVectors[i] = "0" + binaryVectors[i];
        }
      }

      foreach (var binaryVector in binaryVectors) {
        double[] point = new double[dim];
        for (int i = 0; i < binaryVector.Length; i++) {
          if (binaryVector[i] == '0') {
            point[i] = bounds[0];
          } else {
            point[i] = bounds[1];
          }
        }
        hyperCube.Add(point);
      }
      return hyperCube;
    }

    public static double CalculateHypercubeVolume(double[] bounds, int dim) {
      //this is not very smart...
      //return CalculateVolume(GenerateHyperCube(bounds, dim));
      double diff = bounds[1] - bounds[0];

      double volume = diff;
      for (int i = 0; i < dim; i++) {
        volume *= diff;
      }
      return volume;
    }

    public static Fraction CalculateHypercubeVolumeBig(double[] bounds, int dim) {
      double diff = bounds[1] - bounds[0];

      Fraction volume = diff;
      for (int i = 0; i < dim; i++) {
        volume *= diff;
      }
      return volume;
    }
  }
}