1 | function exponent_m = monomialreduction(exponent_m,exponent_p,options,csclasses,LPmodel) |
---|
2 | %MONOMIALREDUCTION Internal function for monomial reduction in SOS programs |
---|
3 | |
---|
4 | % Author Johan Löfberg |
---|
5 | % $Id: monomialreduction.m,v 1.2 2006/09/26 14:28:43 joloef Exp $ |
---|
6 | |
---|
7 | |
---|
8 | % ********************************************** |
---|
9 | % TRIVIAL REDUCTIONS (stupid initial generation) |
---|
10 | % ********************************************** |
---|
11 | mindegrees = min(exponent_p,[],1); |
---|
12 | maxdegrees = max(exponent_p,[],1); |
---|
13 | mindeg = min(sum(exponent_p,2)); |
---|
14 | maxdeg = max(sum(exponent_p,2)); |
---|
15 | if size(exponent_m{1},2)==0 % Stupid case : set(sos(parametric)) |
---|
16 | if options.verbose>0;disp('Initially 1 monomials in R^0');end |
---|
17 | else |
---|
18 | if options.verbose>0;disp(['Initially ' num2str(sum(cellfun('prodofsize',exponent_m)/size(exponent_m{1},2))) ' monomials in R^' num2str(size(exponent_p,2))]);end |
---|
19 | end |
---|
20 | |
---|
21 | for i = 1:length(csclasses) |
---|
22 | t = cputime; |
---|
23 | % THE CODE BELOW IS MESSY TO HANDLE SEVERAL BUGS IN MATLAB |
---|
24 | %too_large_term = any(exponent_m-repmat(maxdegrees/2,size(exponent_m,1),1)>0,2);% DOES NOT HANDLE ODD |
---|
25 | % POLYNIMIALS CORRECTLY |
---|
26 | a1 = full(ceil((1+maxdegrees)/2)); % 6.5.1 in linux freaks on sparse stuff... |
---|
27 | if isempty(a1) |
---|
28 | a1 = zeros(size(maxdegrees)); |
---|
29 | end |
---|
30 | a2 = full(size(exponent_m{i},1)); |
---|
31 | too_large_term = any(exponent_m{i}-repmat(a1,a2,1)>0,2); |
---|
32 | %too_small_term = any(exponent_m-repmat(mindegrees/2,size(exponent_m,1),1)<0,2); |
---|
33 | a1 = full(floor(mindegrees/2)); |
---|
34 | if isempty(a1) |
---|
35 | a1 = zeros(size(mindegrees)); |
---|
36 | end |
---|
37 | a2 = full(size(exponent_m{i},1)); |
---|
38 | too_small_term = any(exponent_m{i}-repmat(a1,a2,1)<0,2);%x^2+xz |
---|
39 | %too_large_sum = any(sum(exponent_m,2)-maxdeg/2>0,2); % DOES NOT HANDLE ODD |
---|
40 | % POLYNIMIALS CORRECTLY |
---|
41 | too_large_sum = any(sum(exponent_m{i},2)-ceil((1+maxdeg)/2)>0,2); |
---|
42 | too_small_sum = any(sum(exponent_m{i},2)-mindeg/2<0,2); |
---|
43 | keep = setdiff1D((1:size(exponent_m{i},1)),find(too_large_term | too_small_term | too_large_sum | too_small_sum)); |
---|
44 | exponent_m{i} = exponent_m{i}(keep,:); |
---|
45 | t = cputime-t; |
---|
46 | end |
---|
47 | if options.verbose>1;disp(['Removing large/small............Keeping ' num2str(sum(cellfun('prodofsize',exponent_m)/size(exponent_m{1},2))) ' monomials (' num2str(t) 'sec)']);end |
---|
48 | |
---|
49 | % ************************************************ |
---|
50 | % Homogenuous? |
---|
51 | % ************************************************ |
---|
52 | if all(sum(exponent_p,2)==sum(exponent_p(1,:))) |
---|
53 | for i = 1:length(csclasses) |
---|
54 | j = csclasses{i}; |
---|
55 | t = cputime; |
---|
56 | exponent_m{i} = exponent_m{i}(sum(exponent_m{i},2)==sum(exponent_p(1,:))/2,:); |
---|
57 | t = cputime-t; |
---|
58 | end |
---|
59 | if options.verbose>1;disp(['Homogenuous form!...............Keeping ' num2str(sum(cellfun('prodofsize',exponent_m)/size(exponent_m{1},2))) ' monomials (' num2str(t) 'sec)']);end |
---|
60 | end |
---|
61 | |
---|
62 | % ************************************************ |
---|
63 | % DIAGONAL CONSISTENCY : MONOMIAL ONLY IN |
---|
64 | % DIAGONAL, CONSTRAINED TO BE ZER0, CAN BE REMOVED |
---|
65 | % ************************************************ |
---|
66 | if (options.sos.inconsistent==1) & ~options.sos.csp |
---|
67 | t = cputime; |
---|
68 | keep = consistent(exponent_m{1},exponent_p); |
---|
69 | exponent_m{1} = exponent_m{1}(keep,:); |
---|
70 | t = cputime-t; |
---|
71 | if options.verbose>0;disp(['Diagonal inconsistensies........Keeping ' num2str(size(exponent_m{1},1)) ' monomials (' num2str(t) 'sec)']);end |
---|
72 | end |
---|
73 | |
---|
74 | % *********************************************************** |
---|
75 | % NEWTON POLYTOPE CHECK |
---|
76 | % *********************************************************** |
---|
77 | if options.sos.newton |
---|
78 | t = cputime; |
---|
79 | [exponent_m,changed,no_lp_solved] = newtonreduce(exponent_m,exponent_p,options,LPmodel); |
---|
80 | t = cputime-t; |
---|
81 | |
---|
82 | if options.verbose>1 | (options.verbose>0 & 1+changed); |
---|
83 | info = ['Newton polytope (' num2str(no_lp_solved) ' LPs)']; |
---|
84 | info = [info repmat('.',1,32-length(info))]; |
---|
85 | disp([info 'Keeping ' num2str(size(exponent_m{end},1)) ' monomials (' num2str(t) 'sec)']); |
---|
86 | end |
---|
87 | end |
---|
88 | |
---|
89 | if (options.sos.inconsistent==2) & ~options.sos.csp |
---|
90 | t = cputime; |
---|
91 | keep = consistent(exponent_m{1},exponent_p); |
---|
92 | exponent_m{1} = exponent_m{1}(keep,:); |
---|
93 | t = cputime-t; |
---|
94 | if options.verbose>1;disp(['Diagonal inconsistensies........Keeping ' num2str(size(exponent_m,1)) ' monomials (' num2str(t) 'sec)']);end |
---|
95 | end |
---|
96 | |
---|
97 | |
---|