source: branches/HeuristicLab.Analysis.AlgorithmBehavior/HeuristicLab.Analysis.AlgorithmBehavior.Analyzers/3.3/LPHull.cs @ 10081

#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 HeuristicLab.Common;

namespace HeuristicLab.Analysis.AlgorithmBehavior.Analyzers {
  /*
   * Calculate Convex Hull using Linear Programming as described in
   * J.B. Rosen et. al., 1989,  Efficient computation of convex hull in Rd
   * and
   * P.M. Pardalos et. al., 1995, 
   * Linear programming approaches to the convex hull problem in Rm 
   * 
   */
  public static class LPHull {
    public static List<double[]> Calculate(double[][] inputs) {
      List<double[]> C = new List<double[]>();
      List<double[]> E = new List<double[]>();

      //Phase 1
      double[][] A = SortA(inputs);

      //Phase 2
      C.Add(A[0]);
      for (int i = 1; i < A.Length; i++) {
        C.Add(A[i]);
        if (!EXT(C, A[i], C.Count - 1)) {
          C.Remove(A[i]);
        }
      }

      //Phase 3
      for (int i = 0; i < C.Count; i++) {
        if (EXT(C, C[i], i)) {
          E.Add(C[i]);
        }
      }

      return E;
    }

    //sort A decreasing by distance to center of polytope
    public static double[][] SortA(double[][] A) {
      int length = A[0].Length;
      double[] maxs = new double[length];
      double[] mins = new double[length];
      double[][] sortedA = new double[A.Length][];

      for (int i = 0; i < maxs.Length; i++) {
        maxs[i] = double.MinValue;
        mins[i] = double.MaxValue;
      }

      for (int i = 0; i < A.Length; i++) {
        for (int j = 0; j < A[i].Length; j++) {
          if (A[i][j] > maxs[j]) maxs[j] = A[i][j];
          if (A[i][j] < mins[j]) mins[j] = A[i][j];
        }
      }

      double[] d = new double[length];
      double[] r = new double[length];

      //calculate center
      for (int i = 0; i < length; i++) {
        d[i] = (maxs[i] + mins[i]) / 2.0;
      }

      //calculate length of sides of rectangle
      for (int i = 0; i < length; i++) {
        r[i] = maxs[i] - mins[i];
      }

      VertexComparer comparer = new VertexComparer();
      comparer.d = d;
      comparer.r = r;
      sortedA = A.OrderByDescending(x => x, comparer).ToArray();

      return sortedA;
    }

    /// <summary>
    /// Checks if alpha is an extreme point (lies on convex hull) of A. 
    /// Returns true if alpha is an extreme point, else false. 
    /// </summary>
    public static bool EXT(List<double[]> A, double[] alpha, int aIdx = -1) {
      alglib.minbleicstate state;
      int N = alpha.Length;
      int K = A.Count;
      double[] init = new double[K];
      double[] lowerBound = new double[K];
      double[] upperBound = new double[K];
      double[] x;
      alglib.minbleicreport rep;
      double[,] c = new double[N + 1, K + 1];
      int[] ct = new int[N + 1]; //init with 0, means equal

      for (int i = 0; i < K; i++) {
        init[i] = 1.0 / K;
        lowerBound[i] = 0.0;
        upperBound[i] = double.MaxValue;
      }

      //last column gets b
      for (int i = 0; i < N; i++) {
        c[i, K] = alpha[i];
      }
      //sum(x) == 1 constraint
      for (int i = 0; i < K + 1; i++) {
        c[N, i] = 1.0;
      }

      //other constraints
      for (int i = 0; i < K; i++) {
        for (int j = 0; j < N; j++) {
          c[j, i] = A[i][j];
        }
      }

      alglib.minbleiccreate(init, out state);
      alglib.minbleicsetbc(state, lowerBound, upperBound);
      alglib.minbleicsetlc(state, c, ct);
      alglib.minbleicsetcond(state, 0.0, 0.0, 0.0, 0);
      alglib.minbleicoptimize(state, Func, null, new Tuple<List<double[]>, double[], int>(A, alpha, aIdx));
      alglib.minbleicresults(state, out x, out rep);

      return x.Count(y => y.IsAlmost(1.0)) == 1;
    }

    private static void Func(double[] x, ref double func, double[] grad, object obj) {
      Tuple<List<double[]>, double[], int> data = obj as Tuple<List<double[]>, double[], int>;
      int aIdx = data.Item3;
      double[] alpha = data.Item2;
      List<double[]> A = data.Item1;

      //well, this seems to work but is probably not correct
      for (int i = 0; i < grad.Length; i++) {
        grad[i] = 0;
      }
      grad[aIdx] = 1;

      func = 0.0;
      for (int i = 0; i < A.Count; i++) {
        for (int j = 0; j < A[i].Length; j++) {
          func += A[i][j] * x[i];
        }
      }
    }
  }
}