#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 .
*/
#endregion
using System;
using System.Collections.Generic;
using System.Linq;
using MIConvexHull;
namespace HeuristicLab.Analysis.AlgorithmBehavior.Analyzers {
public static class ConvexHullMeasures {
//Calculates the volumne of a d-simplex
public static double CalculateSimplexVolume(List 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;
}
//calculates the volume of a convex d-polytope
//decomposition based on delaunay triangulation
public static double CalculateVolume(List 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");
var triangulation = Triangulation.CreateDelaunay(ConvertToVertex(convexHull));
var simplices = ConvertFromTriangulation(triangulation);
foreach (var simplex in simplices) {
volume += CalculateSimplexVolume(simplex.ToList());
}
return volume;
}
private static List ConvertFromTriangulation(ITriangulation> triangulation) {
List results = new List();
foreach (var cell in triangulation.Cells) {
results.Add(cell.Vertices.Select(x => x.Position).ToArray());
}
return results;
}
private static List ConvertToVertex(List data) {
List result = new List();
for (int i = 0; i < data.Count(); i++) {
double[] d = data[i];
DefaultVertex vertex = new DefaultVertex();
vertex.Position = d;
result.Add(vertex);
}
return result;
}
}
}