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