1  #region License Information


2  /* HeuristicLab


3  * Copyright (C) 20022013 Heuristic and Evolutionary Algorithms Laboratory (HEAL)


4  *


5  * This file is part of HeuristicLab.


6  *


7  * HeuristicLab is free software: you can redistribute it and/or modify


8  * it under the terms of the GNU General Public License as published by


9  * the Free Software Foundation, either version 3 of the License, or


10  * (at your option) any later version.


11  *


12  * HeuristicLab is distributed in the hope that it will be useful,


13  * but WITHOUT ANY WARRANTY; without even the implied warranty of


14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the


15  * GNU General Public License for more details.


16  *


17  * You should have received a copy of the GNU General Public License


18  * along with HeuristicLab. If not, see <http://www.gnu.org/licenses/>.


19  */


20  #endregion


21 


22  using System;


23  using System.Collections.Generic;


24  using System.Linq;


25  using HeuristicLab.Common;


26 


27  namespace HeuristicLab.Analysis.AlgorithmBehavior.Analyzers {


28  /*


29  * Calculate Convex Hull using Linear Programming as described in


30  * J.B. Rosen et. al., 1989, Efficient computation of convex hull in Rd


31  * and


32  * P.M. Pardalos et. al., 1995,


33  * Linear programming approaches to the convex hull problem in Rm


34  *


35  */


36  public static class LPHull {


37  public static List<double[]> Calculate(double[][] inputs) {


38  List<double[]> C = new List<double[]>();


39  List<double[]> E = new List<double[]>();


40 


41  //Phase 1


42  double[][] A = SortA(inputs);


43 


44  //Phase 2


45  C.Add(A[0]);


46  for (int i = 1; i < A.Length; i++) {


47  C.Add(A[i]);


48  if (!EXT(C, A[i], C.Count  1)) {


49  C.Remove(A[i]);


50  }


51  }


52 


53  //Phase 3


54  for (int i = 0; i < C.Count; i++) {


55  if (EXT(C, C[i], i)) {


56  E.Add(C[i]);


57  }


58  }


59 


60  return E;


61  }


62 


63  //sort A decreasing by distance to center of polytope


64  public static double[][] SortA(double[][] A) {


65  int length = A[0].Length;


66  double[] maxs = new double[length];


67  double[] mins = new double[length];


68  double[][] sortedA = new double[A.Length][];


69 


70  for (int i = 0; i < maxs.Length; i++) {


71  maxs[i] = double.MinValue;


72  mins[i] = double.MaxValue;


73  }


74 


75  for (int i = 0; i < A.Length; i++) {


76  for (int j = 0; j < A[i].Length; j++) {


77  if (A[i][j] > maxs[j]) maxs[j] = A[i][j];


78  if (A[i][j] < mins[j]) mins[j] = A[i][j];


79  }


80  }


81 


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


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


84 


85  //calculate center


86  for (int i = 0; i < length; i++) {


87  d[i] = (maxs[i] + mins[i]) / 2.0;


88  }


89 


90  //calculate length of sides of rectangle


91  for (int i = 0; i < length; i++) {


92  r[i] = maxs[i]  mins[i];


93  }


94 


95  VertexComparer comparer = new VertexComparer();


96  comparer.d = d;


97  comparer.r = r;


98  sortedA = A.OrderByDescending(x => x, comparer).ToArray();


99 


100  return sortedA;


101  }


102 


103  /// <summary>


104  /// Checks if alpha is an extreme point (lies on convex hull) of A.


105  /// Returns true if alpha is an extreme point, else false.


106  /// </summary>


107  public static bool EXT(List<double[]> A, double[] alpha, int aIdx = 1) {


108  alglib.minbleicstate state;


109  int N = alpha.Length;


110  int K = A.Count;


111  double[] init = new double[K];


112  double[] lowerBound = new double[K];


113  double[] upperBound = new double[K];


114  double[] x;


115  alglib.minbleicreport rep;


116  double[,] c = new double[N + 1, K + 1];


117  int[] ct = new int[N + 1]; //init with 0, means equal


118 


119  for (int i = 0; i < K; i++) {


120  init[i] = 1.0 / K;


121  lowerBound[i] = 0.0;


122  upperBound[i] = double.MaxValue;


123  }


124 


125  //last column gets b


126  for (int i = 0; i < N; i++) {


127  c[i, K] = alpha[i];


128  }


129  //sum(x) == 1 constraint


130  for (int i = 0; i < K + 1; i++) {


131  c[N, i] = 1.0;


132  }


133 


134  //other constraints


135  for (int i = 0; i < K; i++) {


136  for (int j = 0; j < N; j++) {


137  c[j, i] = A[i][j];


138  }


139  }


140 


141  alglib.minbleiccreate(init, out state);


142  alglib.minbleicsetbc(state, lowerBound, upperBound);


143  alglib.minbleicsetlc(state, c, ct);


144  alglib.minbleicsetcond(state, 0.0, 0.0, 0.0, 0);


145  alglib.minbleicoptimize(state, Func, null, new Tuple<List<double[]>, double[], int>(A, alpha, aIdx));


146  alglib.minbleicresults(state, out x, out rep);


147 


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


149  }


150 


151  private static void Func(double[] x, ref double func, double[] grad, object obj) {


152  Tuple<List<double[]>, double[], int> data = obj as Tuple<List<double[]>, double[], int>;


153  int aIdx = data.Item3;


154  double[] alpha = data.Item2;


155  List<double[]> A = data.Item1;


156 


157  //well, this seems to work but is probably not correct


158  for (int i = 0; i < grad.Length; i++) {


159  grad[i] = 0;


160  }


161  grad[aIdx] = 1;


162 


163  func = 0.0;


164  for (int i = 0; i < A.Count; i++) {


165  for (int j = 0; j < A[i].Length; j++) {


166  func += A[i][j] * x[i];


167  }


168  }


169  }


170  }


171  }

