1 | function keep = consistent(exponent_m,exponent_p) |
---|
2 | % CONSISTENT Removes monomials using diagonal inconsistency |
---|
3 | % |
---|
4 | % V = CONSISTENT(V,P) |
---|
5 | % |
---|
6 | % Input |
---|
7 | % P : Scalar SDPVAR object |
---|
8 | % V : Vector with SDPVAR objects |
---|
9 | % |
---|
10 | % Output |
---|
11 | % V : Vector with SDPVAR objects |
---|
12 | % |
---|
13 | % Example: |
---|
14 | % |
---|
15 | % sdpvar x y |
---|
16 | % p = 1+x^4*y^2+x^2*y^4; |
---|
17 | % v = newtonmonoms(p); |
---|
18 | % sdisplay(v) |
---|
19 | % v = consistent(v,p); |
---|
20 | % sdisplay(v) |
---|
21 | % |
---|
22 | % See also NEWTONREDUCE, NEWTONMONOMS, CONGRUENCEBLOCKS |
---|
23 | |
---|
24 | % Author Johan Löfberg |
---|
25 | % $Id: consistent.m,v 1.1 2006/03/30 13:36:20 joloef Exp $ |
---|
26 | |
---|
27 | |
---|
28 | % Default high level call with SDPVAR polynomials |
---|
29 | % Convert to exponent form |
---|
30 | sdpvarout = 0; |
---|
31 | if isa(exponent_m,'sdpvar') |
---|
32 | z = depends(exponent_p); |
---|
33 | z = recover(unique([depends(exponent_p) depends(exponent_m)])); |
---|
34 | [exponent_p,p_base] = getexponentbase(exponent_p,z); |
---|
35 | exponent_m = getexponentbase(exponent_m,z); |
---|
36 | sdpvarout = 1; |
---|
37 | end |
---|
38 | |
---|
39 | n = size(exponent_m,1); |
---|
40 | |
---|
41 | % A bit silly, but I want to keep track |
---|
42 | % on the removed monoms, hence the output |
---|
43 | % from this function is the indicies |
---|
44 | % to the kept monoms |
---|
45 | % |
---|
46 | % Append with indicies |
---|
47 | |
---|
48 | indicies = (1:n)'; |
---|
49 | |
---|
50 | exponent_m = sparse(exponent_m); |
---|
51 | exponent_p = sparse(exponent_p); |
---|
52 | hash = rand(size(exponent_m,2),1); |
---|
53 | exponent_m_hash = exponent_m*hash; |
---|
54 | exponent_2m_hash = (2*exponent_m)*hash; |
---|
55 | exponent_p_hash = (exponent_p)*hash; |
---|
56 | % Check the actual candidates to be removed |
---|
57 | candidates = []; |
---|
58 | for i = 1:n |
---|
59 | m = 2*exponent_m(i,:); |
---|
60 | |
---|
61 | % index_in_p = findrows(exponent_p,m); |
---|
62 | index_in_p = findhash(exponent_p_hash,exponent_2m_hash(i)); |
---|
63 | % if ~isequal(index_in_p,index_in_p2) |
---|
64 | % [index_in_p index_in_p2] |
---|
65 | % end |
---|
66 | if isempty(index_in_p) |
---|
67 | candidates = [candidates;i]; |
---|
68 | end |
---|
69 | end |
---|
70 | |
---|
71 | sums = full(sum(exponent_m,2)); |
---|
72 | for i = 0:max(sums) |
---|
73 | temp{i+1} = exponent_m(sums==i,:); |
---|
74 | end |
---|
75 | maxsum = max(sums); |
---|
76 | |
---|
77 | exponent_2m = 2*exponent_m; |
---|
78 | exponent_m_transpose = exponent_m'; |
---|
79 | if ~isempty(candidates) |
---|
80 | removed = 1; |
---|
81 | while removed |
---|
82 | removelist = []; |
---|
83 | removed=0; |
---|
84 | n = size(exponent_m,1); |
---|
85 | possible = find(ismembc(indicies,candidates)); |
---|
86 | j = 1; |
---|
87 | while j<=length(possible) |
---|
88 | i = possible(j); |
---|
89 | |
---|
90 | % Test this monom |
---|
91 | %m1 = exponent_m(i,:); |
---|
92 | |
---|
93 | %m_squared = 2*m1; |
---|
94 | m_squared = exponent_2m(i,:); |
---|
95 | % Can this be generated else-way, as mj+mk, j,k~=i |
---|
96 | terms = []; |
---|
97 | terms2 = []; |
---|
98 | m_combs = [1:1:i-1 i+1:1:n]; |
---|
99 | m = length(m_combs); |
---|
100 | k = 1; |
---|
101 | |
---|
102 | k = 1; |
---|
103 | %m_squared*hash-exponent_m(m_combs(k),:)*hash |
---|
104 | %full((m_squared*hash- - exponent_m(:,m_combs(k))'*hash) |
---|
105 | while (k<=m) & isempty(terms) |
---|
106 | % m2 = exponent_m(m_combs(k),:); |
---|
107 | m2 = exponent_m_transpose(:,m_combs(k))'; |
---|
108 | x = m_squared-m2; |
---|
109 | % x = 2*m1-m2; |
---|
110 | if all(x>=0) |
---|
111 | terms = find(exponent_m_hash==x*hash); |
---|
112 | end |
---|
113 | k = k+1; |
---|
114 | end |
---|
115 | |
---|
116 | if isempty(terms) |
---|
117 | removelist = [removelist i]; |
---|
118 | removed = 1; |
---|
119 | j = inf; |
---|
120 | end |
---|
121 | j = j+1; |
---|
122 | end |
---|
123 | exponent_m(removelist,:)=[]; |
---|
124 | exponent_m_transpose(:,removelist)=[]; |
---|
125 | exponent_2m(removelist,:)=[]; |
---|
126 | exponent_m_hash(removelist,:)=[]; |
---|
127 | indicies(removelist)=[]; |
---|
128 | end |
---|
129 | end |
---|
130 | keep = indicies; |
---|
131 | |
---|
132 | if sdpvarout |
---|
133 | keep = recovermonoms(exponent_m,z); |
---|
134 | end |
---|
135 | |
---|
136 | return |
---|
137 | |
---|
138 | % **************************** |
---|
139 | % MONOMIAL PRODUCTS |
---|
140 | % **************************** |
---|
141 | % Faster (?) version |
---|
142 | V = symminksum(exponent_m); |
---|
143 | |
---|
144 | |
---|
145 | i = V(:,1)==V(:,2); |
---|
146 | sumlist_diag = V(find(i),:); |
---|
147 | sumlist = V; |
---|
148 | |
---|
149 | % ****************************** |
---|
150 | removed = 1; |
---|
151 | removelist = []; |
---|
152 | slow = 1; |
---|
153 | sumindex = sumlist(:,1:2); |
---|
154 | sumdata = sumlist(:,3:end); |
---|
155 | |
---|
156 | % % Partition data to make searches faster |
---|
157 | merit = sum(sumdata,2) + sum(sumdata>0,2); |
---|
158 | for i = min(merit):max(merit) |
---|
159 | part_sum_data{i+1} = sumdata(find(merit==i),:); |
---|
160 | part_sum_index{i+1} = sumindex(find(merit==i),:); |
---|
161 | end |
---|
162 | |
---|
163 | while removed |
---|
164 | removed = 0; |
---|
165 | removelist = []; |
---|
166 | for i = 1:size(sumlist_diag,1) |
---|
167 | this_monom = sumlist_diag(i,3:end); |
---|
168 | index_in_p = findrows(exponent_p,this_monom); |
---|
169 | if isempty(index_in_p) |
---|
170 | m_this_monom = sum(this_monom) + sum(this_monom>0); |
---|
171 | index_in_m = findrows(part_sum_data{m_this_monom+1},this_monom); |
---|
172 | if ~isempty(index_in_m) & length(index_in_m)==1 |
---|
173 | removed = 1; |
---|
174 | removelist = [removelist sumlist_diag(i,1)]; |
---|
175 | end |
---|
176 | end |
---|
177 | end |
---|
178 | |
---|
179 | for i = 1:length(part_sum_data) |
---|
180 | if ~isempty(part_sum_data{i}) |
---|
181 | index = find(~(ismember(part_sum_index{i}(:,1),removelist) | ismember(part_sum_index{i}(:,2),removelist))); |
---|
182 | part_sum_index{i} = part_sum_index{i}(index,:); |
---|
183 | part_sum_data{i} = part_sum_data{i}(index,:); |
---|
184 | end |
---|
185 | end |
---|
186 | end |
---|
187 | keep = []; |
---|
188 | for i = 1:length(part_sum_data) |
---|
189 | if ~isempty(part_sum_data{i}) |
---|
190 | keep = [keep;part_sum_index{i}]; |
---|
191 | end |
---|
192 | end |
---|
193 | keep = unique(keep); |
---|
194 | %keep = unique(sumindex); |
---|
195 | |
---|
196 | if sdpvarout |
---|
197 | keep = recovermonoms(exponent_m,z); |
---|
198 | end |
---|
199 | |
---|
200 | |
---|
201 | |
---|
202 | |
---|
203 | |
---|
204 | |
---|
205 | function V = symminksum(exponent_m) |
---|
206 | |
---|
207 | n = size(exponent_m,1); |
---|
208 | |
---|
209 | %H1 = kron(speye(n),ones(n,1)); |
---|
210 | %H2 = kron(ones(n,1),speye(n)); |
---|
211 | %V0 = H1*exponent_m + repmat(exponent_m,n,1); |
---|
212 | V0 = kron(exponent_m,ones(n,1)) + repmat(exponent_m,n,1); |
---|
213 | |
---|
214 | ind = []; |
---|
215 | k = 0; |
---|
216 | indicies = []; |
---|
217 | for i = 1:size(exponent_m,1) |
---|
218 | ind = [ind;k+(i:n)'];k = k+n; |
---|
219 | indicies =[indicies;(i:n)' repmat(i,n-i+1,1)]; |
---|
220 | end |
---|
221 | V0 = V0(ind,:); |
---|
222 | V = [indicies V0]; |
---|