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

#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;

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++) {
        if (EXT(C, A[i])) {
          C.Add(A[i]);
        }
      }

      //Phase 3
      for (int i = 0; i < C.Count; i++) {
        //TODO: this is not very efficient
        if (EXT(ExcludeRow(C, i), C[i], i)) {
          E.Add(C[i]);
        }
      }

      return E;
    }

    private static List<double[]> ExcludeRow(List<double[]> matrix, int row) {
      List<double[]> result = new List<double[]>();

      for (int i = 0; i < matrix.Count; i++) {
        if (i != row) {
          var tmpRow = new double[matrix[i].Length];
          result.Add(tmpRow);
          for (int j = 0; j < matrix[i].Length; j++) {
            tmpRow[j] = matrix[i][j];
          }
        }
      }

      return result;
    }

    //sort A decreasing by distance to center of polytope
    public static double[][] SortA(double[][] A) {
      //TODO: sort inplace
      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>
    private static bool EXT(List<double[]> A, double[] alpha, int aIdx = -1) {
      alglib.minbleicstate state;
      double diffstep = 1.0e-6;
      int N = alpha.Length;
      int K = A.Count;


      double[] init = new double[N];
      double[] lowerBound = new double[N];
      double[] upperBound = new double[N];
      double[] x;
      alglib.minbleicreport rep;
      double[,] c = new double[K + 1, N + 1];
      int[] ct = new int[K + 1]; //init with 0, means equal

      for (int i = 0; i < N; i++) {
        //TODO: this should be a feasible solution    
        init[i] = 1.0;
        lowerBound[i] = 0.0;
        upperBound[i] = double.MaxValue;
      }

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

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

      alglib.minbleiccreatef(init, diffstep, 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[]>, int>(A, aIdx));
      alglib.minbleicresults(state, out x, out rep);

      return x.Sum() > 0;
    }

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

      //TODO: I have no idea what I'm doing...
      func = 0.0;
      for (int i = 0; i < A.Count; i++) {
        if (i != aIdx) {
          for (int j = 0; j < A[i].Length; j++) {
            func += A[i][j] * x[i];
          }
        }
      }
    }
  }
}