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

Last change on this file since 10081 was 10081, checked in by ascheibe, 6 years ago

#1886 fixed LPHull, seems to work now

File size: 5.2 KB
Line 
1#region License Information
2/* HeuristicLab
3 * Copyright (C) 2002-2013 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
22using System;
23using System.Collections.Generic;
24using System.Linq;
25using HeuristicLab.Common;
26
27namespace 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}
Note: See TracBrowser for help on using the repository browser.