[37] | 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 |
---|