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