[37] | 1 | function N = congruenceblocks(exponent_m,exponent_p,options,csclasses) |
---|
| 2 | %CONGRUENCEBLOCKS Partitions monomials based on sign symmetry |
---|
| 3 | % |
---|
| 4 | % V = CONGRUENCEBLOCKS(P) |
---|
| 5 | % |
---|
| 6 | % Input |
---|
| 7 | % V : Vector with SDPVAR objects |
---|
| 8 | % P : Scalar SDPVAR object |
---|
| 9 | % |
---|
| 10 | % Output |
---|
| 11 | % V : Cell with SDPVAR objects |
---|
| 12 | % |
---|
| 13 | % Example: |
---|
| 14 | % |
---|
| 15 | % sdpvar x y z |
---|
| 16 | % p = 1+x*y+x^4+y^4+z+z^6; |
---|
| 17 | % v = newtonmonoms(p); |
---|
| 18 | % v = congruenceblocks(v,p); |
---|
| 19 | % sdisplay(v{1}) % Even w.r.t (x,y) |
---|
| 20 | % sdisplay(v{2}) % Odd w.r.t (x,y) |
---|
| 21 | % |
---|
| 22 | % See also NEWTONREDUCE, NEWTONMONOMS, CONSISTENT |
---|
| 23 | |
---|
| 24 | % Author Johan Löfberg |
---|
| 25 | % $Id: congruenceblocks.m,v 1.1 2006/03/30 13:36:20 joloef Exp $ |
---|
| 26 | |
---|
| 27 | sdpvarout = 0; |
---|
| 28 | if isa(exponent_m,'sdpvar') |
---|
| 29 | z = depends(exponent_p); |
---|
| 30 | z = recover(unique([depends(exponent_p) depends(exponent_m)])); |
---|
| 31 | [exponent_p,p_base] = getexponentbase(exponent_p,z); |
---|
| 32 | [m,m_base] = getexponentbase(exponent_m,z); |
---|
| 33 | exponent_m = cell(1);exponent_m{1} = m; |
---|
| 34 | sdpvarout = 1; |
---|
| 35 | end |
---|
| 36 | |
---|
| 37 | if nargin < 3 |
---|
| 38 | options.verbose = 0; |
---|
| 39 | options.sos.congruence = 2; |
---|
| 40 | end |
---|
| 41 | |
---|
| 42 | if nargin < 4 |
---|
| 43 | csclasses = 1; |
---|
| 44 | end |
---|
| 45 | |
---|
| 46 | if ~isempty(exponent_m{1}) & options.sos.congruence>0 & ((size(exponent_p,2)<=16) | options.sos.congruence==1) |
---|
| 47 | |
---|
| 48 | % ********************************************** |
---|
| 49 | % DEFINE CONGRUENCE CLASSES |
---|
| 50 | % ********************************************** |
---|
| 51 | if options.verbose>0;fprintf('Finding symmetries..............');end; |
---|
| 52 | n = size(exponent_p,2); |
---|
| 53 | t = cputime; |
---|
| 54 | switch options.sos.congruence |
---|
| 55 | case 1 |
---|
| 56 | Htemp = eye(n); % CHEAP VERSION; ONLY CHECK IF IT IS EVEN WRT x_i |
---|
| 57 | case 2 |
---|
| 58 | Htemp = dec2decbin(1:2^n-1,n)'; % ALL POSSIBLE COMBINATIONS |
---|
| 59 | otherwise |
---|
| 60 | error('sos.congruence should be 0, 1 or 2') |
---|
| 61 | end |
---|
| 62 | H = Htemp(:,find(~any(rem(exponent_p*Htemp,2),1))); % Find "even" rows |
---|
| 63 | if isempty(H) |
---|
| 64 | N = exponent_m; |
---|
| 65 | if options.verbose>0;disp(['Found no symmetries (' num2str(cputime-t) 'sec)']);end |
---|
| 66 | return |
---|
| 67 | end |
---|
| 68 | |
---|
| 69 | t = cputime-t; |
---|
| 70 | if size(H,2)>=1 |
---|
| 71 | if options.verbose>0 |
---|
| 72 | if size(H,2)>1 |
---|
| 73 | disp(['Found ' num2str(size(H,2)) ' symmetries (' num2str(t) 'sec)']); |
---|
| 74 | else |
---|
| 75 | disp(['Found ' num2str(size(H,2)) ' symmetry (' num2str(t) 'sec)']); |
---|
| 76 | end |
---|
| 77 | end |
---|
| 78 | else |
---|
| 79 | if options.verbose>0;disp(['Found no symmetries (' num2str(t) 'sec)']);end; |
---|
| 80 | end |
---|
| 81 | |
---|
| 82 | % ********************************************** |
---|
| 83 | % CLASSIFY MONMS ACCORDING TO CONGRUENCE CLASSES |
---|
| 84 | % ********************************************** |
---|
| 85 | if size(H,2)>=1 |
---|
| 86 | the_text = 'Partitioning using symmetry.....'; |
---|
| 87 | end |
---|
| 88 | N = cell(0,1); |
---|
| 89 | for cs = 1:length(csclasses) |
---|
| 90 | [ur,j,k]=uniquesafe(mod(exponent_m{cs}*H,2),'rows'); |
---|
| 91 | Ntemp = cell(size(ur,1),1); |
---|
| 92 | temp = []; |
---|
| 93 | for i = 1:length(k) |
---|
| 94 | Ntemp{k(i),1} = [Ntemp{k(i)};exponent_m{cs}(i,:)]; |
---|
| 95 | temp = [temp;size(exponent_m{cs}(i,:),1)]; |
---|
| 96 | end |
---|
| 97 | |
---|
| 98 | for i = 1:length(Ntemp) |
---|
| 99 | N{end+1,1} = Ntemp{i}; |
---|
| 100 | end |
---|
| 101 | end |
---|
| 102 | % ********************************************** |
---|
| 103 | % PRINT SOME RESULTS |
---|
| 104 | % ********************************************** |
---|
| 105 | if size(H,2)>=1 |
---|
| 106 | [uu,ii,oo] = uniquesafe(cellfun('prodofsize',N)/size(N{1},2)); |
---|
| 107 | for i = 1:length(uu) |
---|
| 108 | n_this = length(find(oo==i)); |
---|
| 109 | |
---|
| 110 | the_text = [the_text num2str(uu(i)) 'x' num2str(uu(i)) '(' num2str(n_this) ')' ' ']; |
---|
| 111 | end |
---|
| 112 | end |
---|
| 113 | if options.verbose>0;;disp(the_text);end; |
---|
| 114 | else |
---|
| 115 | N = exponent_m; |
---|
| 116 | end |
---|
| 117 | |
---|
| 118 | if sdpvarout |
---|
| 119 | for i = 1:length(N) |
---|
| 120 | N{i} = recovermonoms(N{i},z); |
---|
| 121 | end |
---|
| 122 | end |
---|