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