1 | function keep = newtonpolytope(exponent_m,exponent_p) |
---|
2 | %NEWTONPOLYTOPE Internal function to remove monimials in SOS programs using Newton polytope |
---|
3 | |
---|
4 | % Author Johan Löfberg |
---|
5 | % $Id: newtonpolytope.m,v 1.1 2006/03/30 13:27:20 joloef Exp $ |
---|
6 | |
---|
7 | |
---|
8 | % WARNING : THIS CODE SUCKS AND IS ONLY USED AS A BACK UP PLAN |
---|
9 | % IF EVERYTHING ELSE FAILS (CRASHING LP SOLVERS ETC) |
---|
10 | |
---|
11 | % ************************************* |
---|
12 | % TRY TO CALCULATE CONVEX HULL |
---|
13 | % ************************************* |
---|
14 | try |
---|
15 | cnvhull = convhulln(full(exponent_p)); |
---|
16 | catch |
---|
17 | keep = 1:size(exponent_m,1); |
---|
18 | return |
---|
19 | end |
---|
20 | |
---|
21 | % *************************************** |
---|
22 | % GET THE UNIQUE POINTS OF IN CONVEX HULL |
---|
23 | % *************************************** |
---|
24 | unique_points = unique(cnvhull); |
---|
25 | |
---|
26 | % *************************************** |
---|
27 | % CALCULATE A POINT IN INTERIOR |
---|
28 | % *************************************** |
---|
29 | p_c = sum(exponent_p(unique_points,:),1)'/length(unique_points); |
---|
30 | |
---|
31 | % *************************************** |
---|
32 | % CALCULATE HYPER-PLANES Ai^T(x-bi)=0 |
---|
33 | % *************************************** |
---|
34 | j = 1; |
---|
35 | A = []; |
---|
36 | for i = 1:size(cnvhull,1) |
---|
37 | X = exponent_p(cnvhull(i,:),:)'; |
---|
38 | y = X(:,1); |
---|
39 | dX = X(:,2:end)-repmat(X(:,1),1,size(X,2)-1); |
---|
40 | Atemp = null(full(dX')); |
---|
41 | % Is this a full-dimensional facet |
---|
42 | if size(Atemp,2)==1 |
---|
43 | direction = (p_c-y)'*Atemp; |
---|
44 | if direction > 0 |
---|
45 | A{j}=-Atemp; |
---|
46 | else |
---|
47 | A{j}=Atemp; |
---|
48 | end |
---|
49 | b{j}=y; |
---|
50 | j=j+1; |
---|
51 | end |
---|
52 | end |
---|
53 | |
---|
54 | % *************************************** |
---|
55 | % CHECK IF MONOMIAL IS IN NEWTON POLYTOPE |
---|
56 | % *************************************** |
---|
57 | if ~isempty(A) |
---|
58 | keep = []; |
---|
59 | for j = 1:size(exponent_m,1) |
---|
60 | inside = 1; |
---|
61 | y = exponent_m(j,:)'; |
---|
62 | if isempty(findrows(exponent_p,y)) % Numerically safe on border |
---|
63 | i = 1; |
---|
64 | while inside & (i<=length(A)) |
---|
65 | inside = inside & ((A{i}'*(y-b{i}))<=1e-9); |
---|
66 | i = i+1; |
---|
67 | end |
---|
68 | end |
---|
69 | if inside |
---|
70 | keep = [keep j]; |
---|
71 | end |
---|
72 | end |
---|
73 | else |
---|
74 | keep = 1:size(exponent_m,1); |
---|
75 | end |
---|