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 |
---|