[37] | 1 | function [base,v] = matrixcoefficients(p,x) |
---|
| 2 | %MATRIXCOEFFICIENTS Extends coefficients, beta |
---|
| 3 | |
---|
| 4 | % Author Johan Löfberg |
---|
| 5 | % $Id: matrixcoefficients.m,v 1.4 2006/08/09 12:14:04 joloef Exp $ |
---|
| 6 | |
---|
| 7 | |
---|
| 8 | % FIX: CURRENTLY UNSTABLE! |
---|
| 9 | |
---|
| 10 | if nargout>1 & (max(size(p))>1) |
---|
| 11 | % error('For matrix inputs, only the coefficients can be returned. Request feature if you need this...'); |
---|
| 12 | end |
---|
| 13 | |
---|
| 14 | % Hack to make sure we get the basis w.r.t all variables and 1 |
---|
| 15 | % This has to be fixed soon (to make robust opt. module fast) |
---|
| 16 | p = p + pi + sum(x)*1e-5; |
---|
| 17 | |
---|
| 18 | if nargin==1 |
---|
| 19 | allvar = depends(p); |
---|
| 20 | xvar = allvar; |
---|
| 21 | x = recover(xvar); |
---|
| 22 | else |
---|
| 23 | xvar = intersect(depends(x),depends(p)); |
---|
| 24 | end |
---|
| 25 | |
---|
| 26 | % Try to debug this! |
---|
| 27 | p = p(:); |
---|
| 28 | base = []; |
---|
| 29 | v = []; |
---|
| 30 | allvar = depends(p);%(ii)); |
---|
| 31 | allvar_recovered = recover(allvar); |
---|
| 32 | t = setdiff(allvar,xvar); |
---|
| 33 | t_recovered = recover(t); |
---|
| 34 | ParametricIndicies = find(ismember(allvar,t)); |
---|
| 35 | map = find(~ismember(allvar,t)); |
---|
| 36 | |
---|
| 37 | for ii = 1:length(p) |
---|
| 38 | pii = p(ii); |
---|
| 39 | [exponent_p,p_base] = getexponentbase(pii,allvar_recovered); |
---|
| 40 | |
---|
| 41 | tempbase = parameterizedbase(pii,[],t_recovered,ParametricIndicies,exponent_p,p_base); |
---|
| 42 | [i,j,k] = unique(full(exponent_p(:,map)),'rows'); |
---|
| 43 | |
---|
| 44 | |
---|
| 45 | V = sparse(1:length(k),k,1,length(tempbase),max(k))'; |
---|
| 46 | |
---|
| 47 | base{ii} = V*tempbase - [pi;repmat(1e-5,length(x),1)]; |
---|
| 48 | |
---|
| 49 | keepthese = j(1:max(k)); |
---|
| 50 | v{ii} = recovermonoms(exponent_p(keepthese,map),x);%recover(xvar)); |
---|
| 51 | end |
---|
| 52 | |
---|
| 53 | |
---|
| 54 | |
---|
| 55 | |
---|
| 56 | |
---|
| 57 | |
---|
| 58 | function p_base_parametric = parameterizedbase(p,z, params,ParametricIndicies,exponent_p,p_base) |
---|
| 59 | |
---|
| 60 | % Check for linear parameterization |
---|
| 61 | parametric_basis = exponent_p(:,ParametricIndicies); |
---|
| 62 | if all(sum(parametric_basis,2)==0) |
---|
| 63 | p_base_parametric = full(p_base(:)); |
---|
| 64 | return |
---|
| 65 | end |
---|
| 66 | if all(sum(parametric_basis,2)<=1) |
---|
| 67 | p_base_parametric = full(p_base(:)); |
---|
| 68 | n = length(p_base_parametric); |
---|
| 69 | ii = []; |
---|
| 70 | vars = []; |
---|
| 71 | js = sum(parametric_basis,1); |
---|
| 72 | for i = 1:size(parametric_basis,2) |
---|
| 73 | if js(i) |
---|
| 74 | j = find(parametric_basis(:,i)); |
---|
| 75 | ii = [ii j(:)']; |
---|
| 76 | vars = [vars repmat(i,1,js(i))]; |
---|
| 77 | end |
---|
| 78 | end |
---|
| 79 | k = setdiff1D(1:n,ii); |
---|
| 80 | if isempty(k) |
---|
| 81 | p_base_parametric = p_base_parametric.*sparse(ii,repmat(1,1,n),params(vars)); |
---|
| 82 | else |
---|
| 83 | pp = params(vars); % Must do this, bug in ML 6.1 (x=sparse(1);x([1 1]) gives different result in 6.1 and 7.0!) |
---|
| 84 | temp = sparse([ii k(:)'],repmat(1,1,n),[pp(:)' ones(1,1,length(k))]); |
---|
| 85 | p_base_parametric = p_base_parametric.*temp; |
---|
| 86 | end |
---|
| 87 | else |
---|
| 88 | % Bummer, nonlinear parameterization sucks... |
---|
| 89 | for i = 1:length(p_base) |
---|
| 90 | j = find(exponent_p(i,ParametricIndicies)); |
---|
| 91 | if ~isempty(j) |
---|
| 92 | temp = p_base(i); |
---|
| 93 | for k = 1:length(j) |
---|
| 94 | if exponent_p(i,ParametricIndicies(j(k)))==1 |
---|
| 95 | temp = temp*params(j(k)); |
---|
| 96 | else |
---|
| 97 | temp = temp*params(j(k))^exponent_p(i,ParametricIndicies(j(k))); |
---|
| 98 | end |
---|
| 99 | end |
---|
| 100 | xx{i} = temp; |
---|
| 101 | else |
---|
| 102 | xx{i} = p_base(i); |
---|
| 103 | end |
---|
| 104 | end |
---|
| 105 | p_base_parametric = stackcell(sdpvar(1,1),xx)'; |
---|
| 106 | end |
---|
| 107 | |
---|
| 108 | |
---|
| 109 | |
---|
| 110 | % |
---|
| 111 | % |
---|
| 112 | % |
---|
| 113 | % |
---|
| 114 | % |
---|
| 115 | % function [base,v] = coefficients(p,x) |
---|
| 116 | % %COEFFICIENTS Extract coefficients and monomials from polynomials |
---|
| 117 | % % |
---|
| 118 | % % [c,v] = COEFFICIENTS(p,x) extracts the coefficents |
---|
| 119 | % % of a polynomial p(x) = c'*v(x) |
---|
| 120 | % % |
---|
| 121 | % % INPUT |
---|
| 122 | % % p : SDPVAR object |
---|
| 123 | % % x : SDPVAR object |
---|
| 124 | % % |
---|
| 125 | % % OUTPUT |
---|
| 126 | % % c : SDPVAR object |
---|
| 127 | % % v : SDPVAR object |
---|
| 128 | % % |
---|
| 129 | % % EXAMPLE |
---|
| 130 | % % sdpvar x y s t |
---|
| 131 | % % p = x^2+x*y*(s+t)+s^2+t^2; % define p(x,y), parameterized with s and t |
---|
| 132 | % % [c,v] = coefficients(p,[x y]); |
---|
| 133 | % % sdisplay([c v]) |
---|
| 134 | % % |
---|
| 135 | % % See also SDPVAR |
---|
| 136 | % |
---|
| 137 | % % Author Johan Löfberg |
---|
| 138 | % % $Id: matrixcoefficients.m,v 1.4 2006/08/09 12:14:04 joloef Exp $ |
---|
| 139 | % |
---|
| 140 | % |
---|
| 141 | % if length(p) > 1%size(p,2) > 1 |
---|
| 142 | % error('Coefficents can only be applied to column vectors'); |
---|
| 143 | % end |
---|
| 144 | % |
---|
| 145 | % allvar = depends(p); |
---|
| 146 | % if nargin==1 |
---|
| 147 | % xvar = allvar; |
---|
| 148 | % x = recover(xvar); |
---|
| 149 | % else |
---|
| 150 | % xvar = depends(x); |
---|
| 151 | % end |
---|
| 152 | % |
---|
| 153 | % pvar = recover(depends(p)); |
---|
| 154 | % |
---|
| 155 | % base = []; |
---|
| 156 | % for i = 1:length(p) |
---|
| 157 | % [bi{i},vi{i}] = coefficientsi(p(i),xvar,pvar,allvar); |
---|
| 158 | % end |
---|
| 159 | % |
---|
| 160 | % % Fix the lengths of the basis to use same basis for all elements |
---|
| 161 | % if length(bi)>1 |
---|
| 162 | % allvars = []; |
---|
| 163 | % for i = 1:length(bi) |
---|
| 164 | % bivar{i} = getvariables(vi{i}); |
---|
| 165 | % if isequal(vi{i}(1),1) |
---|
| 166 | % bivar{i} = [0 bivar{i}]; |
---|
| 167 | % end |
---|
| 168 | % allvars = unique([allvars bivar{i}]); |
---|
| 169 | % end |
---|
| 170 | % v = recover(allvars); |
---|
| 171 | % c = zeros(length(p),length(allvars))'; |
---|
| 172 | % ci = []; |
---|
| 173 | % cj = []; |
---|
| 174 | % cv = []; |
---|
| 175 | % for i = 1:length(bi) |
---|
| 176 | % index = find(ismember(allvars,bivar{i})); |
---|
| 177 | % ci = [ci index]; |
---|
| 178 | % cj = [cj repmat(i,1,length(index))]; |
---|
| 179 | % cv = [cv bi{i}']; |
---|
| 180 | % end |
---|
| 181 | % base = sparse(ci,cj,cv); |
---|
| 182 | % else |
---|
| 183 | % base = bi{1}; |
---|
| 184 | % v = vi{1}; |
---|
| 185 | % end |
---|
| 186 | % |
---|
| 187 | % |
---|
| 188 | % function [base,v] = coefficientsi(p,xvar,pvar,allvar) |
---|
| 189 | % |
---|
| 190 | % % Try to debug this! |
---|
| 191 | % t = setdiff(allvar,xvar); |
---|
| 192 | % [exponent_p,p_base] = getexponentbase(p,pvar); |
---|
| 193 | % ParametricIndicies = find(ismember(allvar,t)); |
---|
| 194 | % % FIX : don't define it here, wait until sparser below. Speed!! |
---|
| 195 | % tempbase = parameterizedbase(p,[],recover(t),ParametricIndicies,exponent_p,p_base); |
---|
| 196 | % [i,j,k] = unique(full(exponent_p(:,find(~ismember(allvar,t)))),'rows'); |
---|
| 197 | % %V = sparse(max(k),length(tempbase)); |
---|
| 198 | % %for i = 1:max(k) |
---|
| 199 | % % V(i,find(k==i)) = 1; |
---|
| 200 | % %end |
---|
| 201 | % V = sparse(1:length(k),k,1,length(tempbase),max(k))'; |
---|
| 202 | % base = V*tempbase; |
---|
| 203 | % if nargout == 2 |
---|
| 204 | % keepthese = j(1:max(k)); |
---|
| 205 | % v = recovermonoms(exponent_p(keepthese,find(~ismember(allvar,t))),recover(xvar)); |
---|
| 206 | % end |
---|
| 207 | % |
---|
| 208 | % |
---|
| 209 | % function p_base_parametric = parameterizedbase(p,z, params,ParametricIndicies,exponent_p,p_base) |
---|
| 210 | % |
---|
| 211 | % % Check for linear parameterization |
---|
| 212 | % parametric_basis = exponent_p(:,ParametricIndicies); |
---|
| 213 | % if all(sum(parametric_basis,2)==0) |
---|
| 214 | % p_base_parametric = full(p_base(:)); |
---|
| 215 | % return |
---|
| 216 | % end |
---|
| 217 | % if all(sum(parametric_basis,2)<=1) |
---|
| 218 | % p_base_parametric = full(p_base(:)); |
---|
| 219 | % n = length(p_base_parametric); |
---|
| 220 | % ii = []; |
---|
| 221 | % vars = []; |
---|
| 222 | % js = sum(parametric_basis,1); |
---|
| 223 | % for i = 1:size(parametric_basis,2) |
---|
| 224 | % if js(i) |
---|
| 225 | % j = find(parametric_basis(:,i)); |
---|
| 226 | % ii = [ii j(:)']; |
---|
| 227 | % vars = [vars repmat(i,1,js(i))]; |
---|
| 228 | % end |
---|
| 229 | % end |
---|
| 230 | % k = setdiff1D(1:n,ii); |
---|
| 231 | % if isempty(k) |
---|
| 232 | % p_base_parametric = p_base_parametric.*sparse(ii,repmat(1,1,n),params(vars)); |
---|
| 233 | % else |
---|
| 234 | % pp = params(vars); % Must do this, bug in ML 6.1 (x=sparse(1);x([1 1]) gives different result in 6.1 and 7.0!) |
---|
| 235 | % p_base_parametric = p_base_parametric.*sparse([ii k(:)'],repmat(1,1,n),[pp(:)' ones(1,1,length(k))]); |
---|
| 236 | % end |
---|
| 237 | % else |
---|
| 238 | % % Bummer, nonlinear parameterization sucks... |
---|
| 239 | % for i = 1:length(p_base) |
---|
| 240 | % j = find(exponent_p(i,ParametricIndicies)); |
---|
| 241 | % if ~isempty(j) |
---|
| 242 | % temp = p_base(i); |
---|
| 243 | % for k = 1:length(j) |
---|
| 244 | % if exponent_p(i,ParametricIndicies(j(k)))==1 |
---|
| 245 | % temp = temp*params(j(k)); |
---|
| 246 | % else |
---|
| 247 | % temp = temp*params(j(k))^exponent_p(i,ParametricIndicies(j(k))); |
---|
| 248 | % end |
---|
| 249 | % end |
---|
| 250 | % xx{i} = temp; |
---|
| 251 | % else |
---|
| 252 | % xx{i} = p_base(i); |
---|
| 253 | % end |
---|
| 254 | % end |
---|
| 255 | % p_base_parametric = stackcell(sdpvar(1,1),xx)'; |
---|
| 256 | % end |
---|