source: proiecte/pmake3d/make3d_original/Make3dSingleImageStanford_version0.1/third_party/opt/yalmip/extras/@lmi/lmi2sedumistruct.m @ 37

Last change on this file since 37 was 37, checked in by (none), 14 years ago

Added original make3d
File size: 13.3 KB
Line 
1function [F_struc,K,KCut] = lmi2sedumistruct(F)
2%lmi2sedumistruct   Internal function: Converts LMI to format needed in SeDuMi
3%
4% % Author Johan Löfberg
5% % $Id: lmi2sedumistruct.m,v 1.21 2006/09/22 08:18:37 joloef Exp $
6
7nvars = yalmip('nvars'); %Needed lot'sa times...
8
9% We first browse to see what we have got and the
10% dimension of F_struc (massive speed improvement)
11%top = 0;
12type_of_constraint = zeros(size(F.clauses,2),1);
13for i = 1:size(F.clauses,2)
14    type_of_constraint(i) = F.clauses{i}.type;
15%    [n,m] = size(F.clauses{i}.data);
16%    [n,m] = size(getbasematrixwithoutcheck(F.clauses{i}.data,0));
17%    top = top+n*m;
18    % Is it a complex linear cone
19 %   if (type_of_constraint(i)==2) & (~isreal(F.clauses{i}.data))
20 %       top = top+n*m; % We will have constraints on Re and Im
21 %   end
22end
23
24F_struc = [];
25
26sdp_con = find(type_of_constraint == 1 | type_of_constraint == 9);
27lin_con = find(type_of_constraint == 2 | type_of_constraint == 12);
28equ_con = find(type_of_constraint == 3);
29qdr_con = find(type_of_constraint == 4);
30rlo_con = find(type_of_constraint == 5);
31
32% SeDuMi struct
33K.f = 0;
34K.l = 0;
35K.q = 0;
36K.r = 0;
37K.s = 0;
38K.rank = [];
39K.dualrank = [];
40K.scomplex = [];
41K.xcomplex = [];
42
43KCut.f = [];
44KCut.l = [];
45KCut.q = [];
46KCut.r = [];
47KCut.s = [];
48
49top = 1;
50localtop = 1;
51% Linear equality constraints
52for i = 1:length(equ_con)
53    constraints = equ_con(i);
54    data = getbase(F.clauses{constraints}.data);
55    [n,m] = size(F.clauses{constraints}.data);
56    % Which variables are needed in this constraint
57    lmi_variables = getvariables(F.clauses{constraints}.data);
58    if isreal(data)
59        ntimesm = n*m; %Just as well pre-calc                   
60    else
61        % Complex constraint, Expand to real and Imag
62        ntimesm = 2*n*m; %Just as well pre-calc
63        data = [real(data);imag(data)];                 
64    end
65    mapX = [1 1+lmi_variables];
66    [ix,jx,sx] = find(data);
67    F_structemp = sparse(ix,mapX(jx),sx,ntimesm,1+nvars);
68    %F_structemp  = spalloc(ntimesm,1+nvars,0);
69    %F_structemp(:,[1 1+lmi_variables(:)'])= data;
70
71    % ...and add them together (efficient for large structures)
72    %   F_struc(top:top+ntimesm-1,:) = F_structemp;     
73    F_struc = [F_struc;F_structemp];
74   
75    if F.clauses{constraints}.cut
76        KCut.f = [KCut.f localtop:localtop+ntimesm-1];
77    end
78   
79    localtop = localtop+ntimesm;           
80    top = top+ntimesm; 
81    K.f = K.f+ntimesm;
82end
83
84% Linear inequality constraints
85localtop = 1;
86F_struc = F_struc';
87% [F_struc,K,KCut] = recursive_lp_fix(F,F_struc,K,KCut,lin_con,nvars,8,1);
88% F_struc = F_struc';
89
90%
91for i = 1:length(lin_con)
92    constraints = lin_con(i);
93    data = getbase(F.clauses{constraints}.data);
94    [n,m] = size(F.clauses{constraints}.data);
95   
96    % Which variables are needed in this constraint
97    lmi_variables = getvariables(F.clauses{constraints}.data);
98   
99    % Convert to real problem
100    if isreal(data)
101        ntimesm = n*m; %Just as well pre-calc                   
102    else
103        % Complex constraint, Expand to real and Imag
104        ntimesm = 2*n*m; %Just as well pre-calc
105        data = [real(data);imag(data)];                 
106    end
107   
108    % Add numerical data to complete problem setup
109    mapX = [1 1+lmi_variables];   
110    [ix,jx,sx] = find(data);   
111    F_structemp = sparse(mapX(jx),ix,sx,1+nvars,ntimesm);
112    F_struc = [F_struc F_structemp];
113   
114    if F.clauses{constraints}.cut
115        KCut.l = [KCut.l localtop:localtop+ntimesm-1];
116    end
117
118    localtop = localtop+ntimesm;       
119    top = top+ntimesm; 
120    K.l = K.l+ntimesm;
121end
122F_struc = F_struc';
123
124[F_struc,K,KCut] = recursive_socp_fix(F,F_struc',K,KCut,qdr_con,nvars,8,1);
125F_struc = F_struc';
126
127
128% Rotated Lorentz cone constraints
129for i = 1:length(rlo_con)
130    constraints = rlo_con(i);
131   
132    [n,m] = size(F.clauses{constraints}.data);
133    ntimesm = n*m; %Just as well pre-calc
134   
135    % Which variables are needed in this constraint
136    lmi_variables = getvariables(F.clauses{constraints}.data);
137   
138    % We allocate the structure blockwise...
139    F_structemp  = spalloc(ntimesm,1+nvars,0);
140    % Add these rows only
141    F_structemp(:,[1 1+lmi_variables(:)'])= getbase(F.clauses{constraints}.data);
142    % ...and add them together (efficient for large structures)
143    %   F_struc(top:top+ntimesm-1,:) = F_structemp;
144    F_struc = [F_struc;F_structemp];
145   
146    top = top+ntimesm;
147    K.r(i) = n;
148end
149
150% Semidefinite  constraints
151% We append the recursively in order to speed up construction
152% of problems with a lot of medium size SDPs
153[F_struc,K,KCut] = recursive_sdp_fix(F,F_struc.',K,KCut,sdp_con,nvars,8,1);
154F_struc = F_struc.';
155
156% Now fix things for the rank constraint
157% This is currently a hack...
158% Should not be in this file
159[rank_variables,dual_rank_variables] = yalmip('rankvariables');
160if ~isempty(rank_variables)   
161    used_in = find(sum(abs(F_struc(:,1+rank_variables)),2));
162    if ~isempty(used_in)
163        if used_in >=1+K.f & used_in < 1+K.l+K.f
164            for i = 1:length(used_in)
165                [ii,jj,kk] = find(F_struc(used_in(i),:));
166                if length(ii)==2 & kk(2)<1
167                    r = floor(kk(1));
168                    var = jj(2)-1;
169                    extstruct = yalmip('extstruct',var);
170                    X = extstruct.arg{1};
171                    if issymmetric(X)
172                        F_structemp = sedumize(X,nvars);
173                    else
174                        error('Only symmetric matrices can be rank constrained.')
175                    end
176                    F_struc = [F_struc;F_structemp];
177                    if isequal(K.s,0)
178                        K.s(1,1) = size(extstruct.arg{1},1);
179                    else
180                        K.s(1,end+1) = size(extstruct.arg{1},1);
181                    end
182                    K.rank(1,end+1) = min(r,K.s(end));
183                else
184                    error('This rank constraint is not supported (only supports rank(X) < r)')
185                end
186            end
187            % Remove the nonlinear operator constraints
188           
189            F_struc(used_in,:) = [];
190            K.l = K.l - length(used_in);
191        else
192            error('You have added a rank constraint on an equality constraint, or a scalar expression?!')
193        end
194    end
195end
196if ~isempty(dual_rank_variables)
197    used_in = find(sum(abs(F_struc(:,1+dual_rank_variables)),2));
198    if ~isempty(used_in)
199        if used_in >=1+K.f & used_in < 1+K.l+K.f
200            for i = 1:length(used_in)
201                [ii,jj,kk] = find(F_struc(used_in(i),:));
202                if length(ii)==2 & kk(2)<1
203                    r = floor(kk(1));
204                    var = jj(2)-1;
205                    extstruct = yalmip('extstruct',var);
206                    X = extstruct.arg{1};
207                    id = getlmiid(X);
208                    inlist=getlmiid(F);
209                    index=find(id==inlist);
210                    if ~isempty(index)                                           
211                        K.rank(1,index) = min(r,K.s(index));
212                    end
213                else
214                    error('This rank constraint is not supported (only supports rank(X) < r)')
215                end
216            end
217            % Remove the nonlinear operator constraints
218            F_struc(used_in,:) = [];
219            K.l = K.l - length(used_in);
220        else
221            error('You have added a rank constraint on an equality constraint, or a scalar expression?!')
222        end
223    end
224end
225
226function F_structemp = sedumize(Fi,nvars)
227Fibase = getbase(Fi);
228[n,m] = size(Fi);
229ntimesm = n*m;
230lmi_variables = getvariables(Fi);
231[ix,jx,sx] = find(Fibase);
232mapX = [1 1+lmi_variables];
233F_structemp = sparse(ix,mapX(jx),sx,ntimesm,1+nvars);
234
235function [F_struc,K,KCut] = recursive_lp_fix(F,F_struc,K,KCut,lp_con,nvars,maxnlp,startindex)
236
237% Check if we should recurse
238if length(lp_con)>2*maxnlp
239    % recursing costs, so do 4 in one step
240    ind = 1+ceil(length(lp_con)*(0:0.25:1));
241    [F_struc1,K,KCut] = recursive_lp_fix(F,[],K,KCut,lp_con(ind(1):ind(2)-1),nvars,maxnlp,startindex+ind(1)-1);
242    [F_struc2,K,KCut] = recursive_lp_fix(F,[],K,KCut,lp_con(ind(2):ind(3)-1),nvars,maxnlp,startindex+ind(2)-1);
243    [F_struc3,K,KCut] = recursive_lp_fix(F,[],K,KCut,lp_con(ind(3):ind(4)-1),nvars,maxnlp,startindex+ind(3)-1);
244    [F_struc4,K,KCut] = recursive_lp_fix(F,[],K,KCut,lp_con(ind(4):ind(5)-1),nvars,maxnlp,startindex+ind(4)-1);
245    F_struc = [F_struc F_struc1 F_struc2 F_struc3 F_struc4];
246    return
247elseif length(lp_con)>maxnlp
248    mid = ceil(length(lp_con)/2);
249    [F_struc1,K,KCut] = recursive_lp_fix(F,[],K,KCut,lp_con(1:mid),nvars,maxnlp,startindex);
250    [F_struc2,K,KCut] = recursive_lp_fix(F,[],K,KCut,lp_con(mid+1:end),nvars,maxnlp,startindex+mid);
251    F_struc = [F_struc F_struc1 F_struc2];
252    return
253end
254
255oldF_struc = F_struc;
256F_struc = [];
257for i = 1:length(lp_con)
258    constraints = lp_con(i);
259    Fi = F.clauses{constraints}.data;
260    Fibase = getbase(Fi);
261    [n,m] = size(Fi);
262   
263    %ntimesm = n*m; %Just as well pre-calc
264   
265   
266    % Convert to real problem
267    if isreal(Fibase)
268        ntimesm = n*m; %Just as well pre-calc                   
269    else
270        % Complex constraint, Expand to real and Imag
271        ntimesm = 2*n*m; %Just as well pre-calc
272        Fibase = [real(Fibase);imag(Fibase)];                   
273    end
274
275    % Which variables are needed in this constraint
276    lmi_variables = getvariables(Fi);
277    mapX = [1 1+lmi_variables];
278   
279    [ix,jx,sx] = find(Fibase);
280   
281    F_structemp = sparse(mapX(jx),ix,sx,1+nvars,ntimesm);
282    F_struc = [F_struc F_structemp];
283       
284    if F.clauses{constraints}.cut
285        KCut.l = [KCut.l i+startindex-1];
286    end
287   
288    K.l(i+startindex-1) = n;   
289end
290K.l = sum(K.l);
291F_struc = [oldF_struc F_struc];
292
293
294function [F_struc,K,KCut] = recursive_sdp_fix(F,F_struc,K,KCut,sdp_con,nvars,maxnsdp,startindex)
295
296% Check if we should recurse
297if length(sdp_con)>2*maxnsdp
298    % recursing costs, so do 4 in one step
299    ind = 1+ceil(length(sdp_con)*(0:0.25:1));
300    [F_struc1,K,KCut] = recursive_sdp_fix(F,[],K,KCut,sdp_con(ind(1):ind(2)-1),nvars,maxnsdp,startindex+ind(1)-1);
301    [F_struc2,K,KCut] = recursive_sdp_fix(F,[],K,KCut,sdp_con(ind(2):ind(3)-1),nvars,maxnsdp,startindex+ind(2)-1);
302    [F_struc3,K,KCut] = recursive_sdp_fix(F,[],K,KCut,sdp_con(ind(3):ind(4)-1),nvars,maxnsdp,startindex+ind(3)-1);
303    [F_struc4,K,KCut] = recursive_sdp_fix(F,[],K,KCut,sdp_con(ind(4):ind(5)-1),nvars,maxnsdp,startindex+ind(4)-1);
304    F_struc = [F_struc F_struc1 F_struc2 F_struc3 F_struc4];
305    return
306elseif length(sdp_con)>maxnsdp
307    mid = ceil(length(sdp_con)/2);
308    [F_struc1,K,KCut] = recursive_sdp_fix(F,[],K,KCut,sdp_con(1:mid),nvars,maxnsdp,startindex);
309    [F_struc2,K,KCut] = recursive_sdp_fix(F,[],K,KCut,sdp_con(mid+1:end),nvars,maxnsdp,startindex+mid);
310    F_struc = [F_struc F_struc1 F_struc2];
311    return
312end
313
314oldF_struc = F_struc;
315F_struc = [];
316for i = 1:length(sdp_con)
317    constraints = sdp_con(i);
318    Fi = F.clauses{constraints}.data;
319    Fibase = getbase(Fi);
320    [n,m] = size(Fi);
321    ntimesm = n*m; %Just as well pre-calc
322
323    % Which variables are needed in this constraint
324    lmi_variables = getvariables(Fi);
325    mapX = [1 1+lmi_variables];
326   
327    [ix,jx,sx] = find(Fibase);
328   
329    F_structemp = sparse(mapX(jx),ix,sx,1+nvars,ntimesm);
330    F_struc = [F_struc F_structemp];
331
332    if F.clauses{constraints}.cut
333        KCut.s = [KCut.s i+startindex-1];
334    end
335    K.s(i+startindex-1) = n;
336    K.rank(i+startindex-1) = n;
337    K.dualrank(i+startindex-1) = n;
338    % Check for a complex structure
339    if ~isreal(F_structemp)
340        K.scomplex = [K.scomplex i+startindex-1];
341    end
342end
343F_struc = [oldF_struc F_struc];
344
345
346
347
348
349function [F_struc,K,KCut] = recursive_socp_fix(F,F_struc,K,KCut,qdr_con,nvars,maxnsocp,startindex);
350
351% Check if we should recurse
352if length(qdr_con)>2*maxnsocp
353    % recursing costs, so do 4 in one step
354    ind = 1+ceil(length(qdr_con)*(0:0.25:1));
355    [F_struc1,K,KCut] = recursive_socp_fix(F,[],K,KCut,qdr_con(ind(1):ind(2)-1),nvars,maxnsocp,startindex+ind(1)-1);
356    [F_struc2,K,KCut] = recursive_socp_fix(F,[],K,KCut,qdr_con(ind(2):ind(3)-1),nvars,maxnsocp,startindex+ind(2)-1);
357    [F_struc3,K,KCut] = recursive_socp_fix(F,[],K,KCut,qdr_con(ind(3):ind(4)-1),nvars,maxnsocp,startindex+ind(3)-1);
358    [F_struc4,K,KCut] = recursive_socp_fix(F,[],K,KCut,qdr_con(ind(4):ind(5)-1),nvars,maxnsocp,startindex+ind(4)-1);
359    F_struc = [F_struc F_struc1 F_struc2 F_struc3 F_struc4];
360    return
361elseif length(qdr_con)>maxnsocp
362    mid = ceil(length(qdr_con)/2);
363    [F_struc1,K,KCut] = recursive_socp_fix(F,[],K,KCut,qdr_con(1:mid),nvars,maxnsocp,startindex);
364    [F_struc2,K,KCut] = recursive_socp_fix(F,[],K,KCut,qdr_con(mid+1:end),nvars,maxnsocp,startindex+mid);
365    F_struc = [F_struc  F_struc1  F_struc2];
366    return
367end
368
369% second order cone constraints
370for i = 1:length(qdr_con)
371    constraints = qdr_con(i);
372
373    [n,m] = size(F.clauses{constraints}.data);
374    ntimesm = n*m; %Just as well pre-calc
375
376    % Which variables are needed in this constraint
377    lmi_variables = getvariables(F.clauses{constraints}.data);
378
379    data = getbase(F.clauses{constraints}.data);
380    if isreal(data)
381        mapX = [1 1+lmi_variables];
382        [ix,jx,sx] = find(data);
383        F_structemp = sparse(mapX(jx),ix,sx,1+nvars,ntimesm);
384    else
385        n = n+(n-1);
386        ntimesm = n*m;
387        F_structemp  = spalloc(ntimesm,1+nvars,0);
388        data = [data(1,:);real(data(2:end,:));imag(data(2:end,:))];
389        F_structemp(:,[1 1+lmi_variables(:)'])= data;
390        F_structemp = F_structemp';
391    end
392    % ...and add them together (efficient for large structures)
393    F_struc = [F_struc F_structemp];
394    K.q(i+startindex-1) = n;
395end
396
397
398
399
400
401
Note: See TracBrowser for help on using the repository browser.