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

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

Added original make3d
File size: 13.2 KB
Line 
1function [problem,integer_variables,binary_variables,parametric_variables,uncertain_variables,quad_info] = categorizeproblem(F,P,h,relax,parametric,evaluation)
2%categorizeproblem          Internal function: tries to determine the type of optimization problem
3
4% Author Johan Löfberg
5% $Id: categorizeproblem.m,v 1.12 2006/09/22 11:57:38 joloef Exp $
6
7Counter = size(F.clauses,2);
8Ftype = zeros(Counter,1);
9
10real_data = 1;
11int_data  = 0;
12bin_data  = 0;
13par_data  = 0;
14
15poly_constraint = 0;
16bilin_constraint = 0;
17sigm_constraint = 0;
18rank_constraint = 0;
19rank_objective = 0;
20
21integer_variables = [];
22binary_variables = [];
23parametric_variables = [];
24kyp_prob  = 0;
25
26% ***********************************************
27% Setup an empty problem definition
28% ***********************************************
29problem.objective.linear = 0;
30problem.objective.quadratic.convex = 0;
31problem.objective.quadratic.nonconvex = 0;
32problem.objective.polynomial = 0;
33problem.objective.maxdet = 0;
34problem.objective.sigmonial = 0;
35
36problem.constraint.equalities.linear     = 0;
37problem.constraint.equalities.quadratic  = 0;
38problem.constraint.equalities.polynomial = 0;
39problem.constraint.equalities.sigmonial  = 0;
40
41problem.constraint.inequalities.elementwise.linear = 0;
42problem.constraint.inequalities.elementwise.quadratic.convex = 0;
43problem.constraint.inequalities.elementwise.quadratic.nonconvex = 0;
44problem.constraint.inequalities.elementwise.sigmonial = 0;
45problem.constraint.inequalities.elementwise.polynomial = 0;
46
47problem.constraint.inequalities.semidefinite.linear = 0;
48problem.constraint.inequalities.semidefinite.quadratic = 0;
49problem.constraint.inequalities.semidefinite.polynomial = 0;
50problem.constraint.inequalities.semidefinite.sigmonial = 0;
51
52problem.constraint.inequalities.rank = 0;
53
54problem.constraint.inequalities.secondordercone = [];
55problem.constraint.inequalities.rotatedsecondordercone = [];
56
57problem.constraint.integer = 0;
58problem.constraint.binary = 0;
59
60problem.complex = 0;
61problem.parametric = parametric;
62problem.evaluation = evaluation;
63
64% ********************************************************
65% Make a list of all globally available discrete variables
66% ********************************************************
67integer_variables   = yalmip('intvariables');
68binary_variables    = yalmip('binvariables');
69uncertain_variables = yalmip('uncvariables');
70for i = 1:Counter
71    switch F.clauses{i}.type
72        case 7
73            integer_variables = union(integer_variables,getvariables(F.clauses{i}.data));
74        case 8
75            binary_variables = union(binary_variables,getvariables(F.clauses{i}.data));
76        case 13
77            parametric_variables = union(parametric_variables,getvariables(F.clauses{i}.data));
78        otherwise
79    end
80end
81
82% ********************************************************
83% Logarithmic objective?
84% ********************************************************
85problem.objective.maxdet = ~isempty(P);
86
87% ********************************************************
88% Rank variables
89% ********************************************************
90rank_variables = yalmip('rankvariables');
91any_rank_variables = ~isempty(rank_variables);
92
93% ********************************************************
94% Make a list of all globally available nonlinear variables
95% ********************************************************
96[monomtable,variabletype] = yalmip('monomtable');
97
98linear_variables = find(variabletype==0);
99nonlinear_variables = find(variabletype~=0);
100sigmonial_variables = find(variabletype==4);
101
102allvars = getvariables(F);
103any_nonlinear_variables =~isempty(find(ismembc(nonlinear_variables,allvars)));
104any_discrete_variables = ~isempty(integer_variables) | ~isempty(binary_variables);
105
106for i = 1:Counter
107   
108    Fi = F.clauses{i};
109    % Only real-valued data?
110    real_data = real_data & isreal(F.clauses{i}.data);
111   
112    % Any discrete variables used
113    if any_discrete_variables
114        Fvar = getvariables(Fi.data);
115        int_data = int_data | any(ismembc(Fvar,integer_variables));
116        bin_data = bin_data | any(ismembc(Fvar,binary_variables));
117        par_data = par_data | any(ismembc(Fvar,parametric_variables));
118    end
119   
120    if any_rank_variables       
121        rank_constraint = rank_constraint | any(ismember(getvariables(Fi.data),rank_variables));
122    end 
123   
124    if ~any_nonlinear_variables % No nonlinearly parameterized constraints
125       
126        switch Fi.type
127        case {1,9}
128            problem.constraint.inequalities.semidefinite.linear = 1;
129        case 2
130            problem.constraint.inequalities.elementwise.linear = 1;
131        case 3
132            problem.constraint.equalities.linear = 1;
133        case 4
134            problem.constraint.inequalities.secondordercone = 1;
135        case 5
136            problem.constraint.inequalities.rotatedsecondordercone = 1;           
137        otherwise
138        end
139    else
140        % Can be nonlinear stuff
141        vars = getvariables(Fi.data);
142        usednonlins = find(ismembc(nonlinear_variables,vars));
143        if ~isempty(usednonlins)
144            usedsigmonials = find(ismember(sigmonial_variables,vars));
145            if ~isempty(usedsigmonials)
146                switch Fi.type
147                case 1
148                    problem.constraint.inequalities.semidefinite.sigmonial = 1;
149                case 2
150                    problem.constraint.inequalities.elementwise.sigmonial = 1;
151                case 3
152                    problem.constraint.equalities.sigmonial = 1;
153                case 4
154                    error('Sigmonial SOCP not supported');                     
155                case 5
156                    error('Sigmonial RSOCP not supported');                   
157                otherwise
158                    error('Report bug in problem classification (sigmonial constraint)');
159                end
160            else
161                deg = degree(Fi.data);
162                switch deg
163                   
164                case 1
165                    switch Fi.type
166                    case 1
167                        problem.constraint.inequalities.semidefinite.linear = 1;
168                    case 2
169                        problem.constraint.inequalities.elementwise.linear = 1;
170                    case 3
171                        problem.constraint.equalities.linear = 1;
172                    case 4
173                        problem.constraint.inequalities.secondordercone = 1;
174                    case 5
175                        problem.constraint.inequalities.rotatedsecondordercone = 1;                               
176                    otherwise
177                        error('Report bug in problem classification (linear constraint)');
178                    end
179                case 2
180                    switch Fi.type
181                    case 1
182                        problem.constraint.inequalities.semidefinite.quadratic = 1;
183                    case 2
184                        % FIX : This should be re-used from
185                        % convertconvexquad
186                        convex = 1;
187                        f = Fi.data;f = f(:);
188                        ii = 1;
189                        while convex & ii<=length(f)
190                            [Q,c,f,x,info] = quaddecomp(f(ii));
191
192                            if info | any(eig(Q) > 0)
193                                convex = 0;
194                            end
195                            ii= ii + 1;
196                        end
197                        if convex
198                            problem.constraint.inequalities.elementwise.quadratic.convex = 1;
199                        else
200                            problem.constraint.inequalities.elementwise.quadratic.nonconvex = 1;
201                        end
202                    case 3
203                        problem.constraint.equalities.quadratic = 1;
204                    case 4
205                        error                               
206                    case 5
207                        error                               
208                    otherwise
209                        error('Report bug in problem classification (quadratic constraint)');
210                    end
211                otherwise
212                    switch Fi.type
213                    case 1
214                        problem.constraint.inequalities.semidefinite.polynomial = 1;
215                    case 2
216                        problem.constraint.inequalities.elementwise.polynomial = 1;
217                    case 3
218                        problem.constraint.equalities.polynomial = 1;
219                    case 4
220                        %   problem.constraint.inequalities.secondordercone = 1;
221                    case 5
222                        %   problem.constraint.inequalities.rotatedsecondordercone = 1;
223                    otherwise
224                        error('Report bug in problem classification (polynomial constraint)');
225                    end
226                   
227                end
228            end
229        else
230            switch Fi.type
231            case 1
232                problem.constraint.inequalities.semidefinite.linear = 1;
233            case 2
234                problem.constraint.inequalities.elementwise.linear = 1;
235            case 3
236                problem.constraint.equalities.linear = 1;
237            case 4
238                problem.constraint.inequalities.secondordercone = 1;                   
239            case 5
240                problem.constraint.inequalities.rotatedsecondordercone = 1;                                                   
241            case 7
242                problem.constraint.integer = 1;
243            case 8
244                problem.constraint.binary = 1;                     
245            otherwise
246                error('Report bug in problem classification (linear constraint)');
247            end
248        end
249    end
250end
251
252if int_data   
253    problem.constraint.integer = 1;
254end
255if bin_data
256    problem.constraint.binary = 1;
257end
258if ~real_data
259    problem.complex = 1;
260end
261if rank_constraint
262    problem.constraint.inequalities.rank = 1;
263end
264if ~isempty(uncertain_variables)
265    problem.uncertain = 1;
266end
267
268if (relax==1) | (relax==2)
269    problem.constraint.integer = 0;
270    problem.constraint.binary = 0;       
271    int_data  = 0;
272    bin_data  = 0;
273%    integer_variables = [];
274%    binary_variables = [];
275end
276if (relax==1) | (relax==3)   
277    problem.constraint.equalities.linear = problem.constraint.equalities.linear | problem.constraint.equalities.quadratic | problem.constraint.equalities.polynomial | problem.constraint.equalities.sigmonial;
278    problem.constraint.equalities.quadratic = 0;
279    problem.constraint.equalities.polynomial = 0;
280    problem.constraint.equalities.sigmonial = 0;
281   
282    problem.constraint.inequalities.elementwise.linear = problem.constraint.inequalities.elementwise.linear | problem.constraint.inequalities.elementwise.quadratic.convex | problem.constraint.inequalities.elementwise.quadratic.nonconvex | problem.constraint.inequalities.elementwise.sigmonial | problem.constraint.inequalities.elementwise.polynomial;
283    problem.constraint.inequalities.elementwise.quadratic.convex = 0;
284    problem.constraint.inequalities.elementwise.quadratic.nonconvex = 0;
285    problem.constraint.inequalities.elementwise.sigmonial   = 0;
286    problem.constraint.inequalities.elementwise.polynomial  = 0;   
287   
288    problem.constraint.inequalities.semidefinite.linear =  problem.constraint.inequalities.semidefinite.linear | problem.constraint.inequalities.semidefinite.quadratic | problem.constraint.inequalities.semidefinite.polynomial | problem.constraint.inequalities.semidefinite.sigmonial;   
289    problem.constraint.inequalities.semidefinite.quadratic  = 0;
290    problem.constraint.inequalities.semidefinite.polynomial = 0;
291    problem.constraint.inequalities.semidefinite.sigmonial  = 0;
292   
293    poly_constraint = 0;
294    bilin_constraint = 0;
295    sigm_constraint = 0;   
296end
297
298
299% Analyse the objective function
300quad_info = [];
301if (~isempty(h)) & ~is(h,'linear') &~(relax==1) &~(relax==3)
302    if any(ismember(getvariables(h),sigmonial_variables))
303        problem.objective.sigmonial = 1;
304    else
305        [Q,c,f,x,info] = quaddecomp(h);
306        if ~isreal(Q) % Numerical noise common on imaginary parts
307            Qr = real(Q);
308            Qi = imag(Q);
309            Qr(abs(Qr)<1e-10) = 0;
310            Qi(abs(Qi)<1e-10) = 0;
311            cr = real(c);
312            ci = imag(c);
313            cr(abs(cr)<1e-10) = 0;
314            ci(abs(ci)<1e-10) = 0;
315            Q = Qr + sqrt(-1)*Qi;
316            c = cr + sqrt(-1)*ci;           
317        end
318        if info==0         
319            [R,p]=chol(Q);
320            if p~=0
321                Q = full(Q);
322                if min(eig(Q))>=-1e-10
323                    p=0;     
324                    try
325                        [U,S,V]=svd(Q);
326                    catch
327                        [U,S,V]=svd(full(Q));
328                    end
329                    i = find(diag(S)>1e-10);
330                    R = sqrt(S(1:max(i),1:max(i)))*V(:,1:max(i))';
331                end
332            end
333            if p==0
334                problem.objective.quadratic.convex = 1;
335            else
336                problem.objective.quadratic.nonconvex = 1;
337            end
338            quad_info.Q = Q;
339            quad_info.c = c;
340            quad_info.f = f;
341            quad_info.x = x;
342            quad_info.R = R;
343            quad_info.p = p;
344        else
345            problem.objective.polynomial = 1;
346        end
347    end
348else   
349    problem.objective.linear = ~isempty(h);
350end
Note: See TracBrowser for help on using the repository browser.