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

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

Added original make3d
File size: 8.4 KB
Line 
1function [F_expand,failure,cause] = expandrecursive(variable,F_expand,extendedvariables,monomtable,where,level,options,method,extstruct,goal_vexity)
2global DUDE_ITS_A_GP ALREADY_MODELLED REMOVE_THESE_IN_THE_END MARKER_VARIABLES OPERATOR_IN_POLYNOM
3cause = '';
4failure = 0;
5
6if  ~alreadydone(getvariables(variable))
7    % Request epigraph (or hypograph) for this variable
8    [F_graph,properties,arguments,fcn] = model(variable,method,options,extstruct);
9
10    % This is useful in MPT
11    if ~isempty(F_graph)
12       F_graph = tag(F_graph,['Expansion of ' fcn]);
13    end
14   
15    % Bit of a hack (CPLEX sucks on some trivial problems, so we have to track
16    % these simple models and fix them in the end)
17    if ~isempty(properties)
18        if isfield(properties,'extra')
19            if strcmp(properties.extra,'marker')
20                MARKER_VARIABLES = [MARKER_VARIABLES getvariables(variable)];
21            end
22        end
23    end
24
25    % Now check that this operator actually is convex/concave/milp
26    if isequal(method,'graph')
27
28        switch properties.convexity
29            case {'convex','concave'}
30                failure = ~isequal(properties.convexity,goal_vexity);
31            case 'milp'
32                if options.allowmilp
33                    failure = 0;
34                    method = 'milp';
35                else
36                    failure = 1;
37                end
38            case {'none','exact'}
39                failure = 0;
40        end
41
42        if failure & options.allowmilp
43            [F_graph,properties,arguments] = model(variable,'milp',options);           
44            % This is useful in MPT
45            if ~isempty(F_graph)
46               F_graph = tag(F_graph,['Expansion of ' fcn]);
47            end
48            if isempty(F_graph)
49                cause = [cause ', MILP model not available.'];
50                return
51            else
52                failure = 0;
53                method = 'milp';
54            end
55        elseif failure
56            cause = ['Expected ' 'goal_vexity' ' function in '  where ' at level ' num2str(level)];
57        end
58    else
59        if isempty(F_graph)
60            cause = ['MILP model not available in ' where ' at level ' num2str(level)];
61            failure = 1;
62            return
63        else
64            failure = 0;
65            method = 'milp';
66        end
67    end
68
69    % We save variable models for future use
70    if failure == 0
71        done = save_model_expansion(method,F_graph,properties);
72        if done
73            return
74        end
75    end
76
77    % Now we might have to recurse
78    arguments = arguments(:);
79
80    [ix,jx,kx] = find(monomtable(getvariables(variable),:));
81    if ~isempty(jx) % Bug in 6.1
82        if any(kx>1)
83            OPERATOR_IN_POLYNOM = [OPERATOR_IN_POLYNOM extendedvariables(jx(find(kx>1)))];
84        end
85    end
86
87    % Small pre-processing to speed-up large-scale problems (subsref sloooow)
88    % with only linear arguments (such as norm(Ax-b) problems)
89    if isa(arguments,'sdpvar')
90        do_not_check_nonlinearity = is(arguments,'linear');
91        if do_not_check_nonlinearity
92            allvariables = getvariables(arguments);
93            fullbasis = getbase(arguments);
94            fullbasis = fullbasis(:,2:end);
95            fullbasis_transpose = fullbasis';
96        end
97    else
98        do_not_check_nonlinearity = 0;
99    end
100
101    j = 1;
102    % Ok, here goes the actual recursive code 
103    while j<=length(arguments) & ~failure
104
105        if do_not_check_nonlinearity
106%            usedvariables = find(fullbasis(j,:));
107            usedvariables = find(fullbasis_transpose(:,j));
108            expressionvariables = allvariables(usedvariables);
109        else
110            expression = arguments(j);
111            expressionvariables = unique([depends(expression) getvariables(expression)]);
112        end
113        index_in_expression = find(ismembc(expressionvariables,extendedvariables));
114
115        if ~isempty(index_in_expression)
116            for i = index_in_expression
117                if ~alreadydone(expressionvariables(i))
118                    if do_not_check_nonlinearity
119                      %  basis = fullbasis(j,expressionvariables(index_in_expression(i)));
120                       % basis = fullbasis(j,usedvariables(index_in_expression(i)));
121                        basis = fullbasis(j,usedvariables((i)));
122                    else
123                        basis = getbasematrix(expression,expressionvariables(i));
124                    end
125
126                    go_convex1 = (basis > 0) &  isequal(goal_vexity,'convex') & isequal(properties.monotonicity,'increasing');
127                    go_convex2 = (basis <= 0) &  isequal(goal_vexity,'convex') & isequal(properties.monotonicity,'decreasing');
128                    go_convex3 = (basis <= 0) &  isequal(goal_vexity,'concave') & isequal(properties.monotonicity,'increasing');
129                    go_convex4 = (basis > 0) &  isequal(goal_vexity,'concave') & isequal(properties.monotonicity,'decreasing');
130
131                    go_concave1 = (basis > 0) &  isequal(goal_vexity,'convex') & isequal(properties.monotonicity,'decreasing');
132                    go_concave2 = (basis <= 0) &  isequal(goal_vexity,'convex') & isequal(properties.monotonicity,'increasing');
133                    go_concave3 = (basis <= 0) &  isequal(goal_vexity,'concave') & isequal(properties.monotonicity,'decreasing');
134                    go_concave4 = (basis > 0) &  isequal(goal_vexity,'concave') & isequal(properties.monotonicity,'increasing');
135
136                    if go_convex1 | go_convex2 | go_convex3 | go_convex4
137                        [F_expand,failure,cause] = expandrecursive(recover(expressionvariables(i)),F_expand,extendedvariables,monomtable,where,level+1,options,method,[],'convex');
138                    elseif go_concave1 | go_concave2 | go_concave3 | go_concave4
139                        [F_expand,failure,cause] = expandrecursive(recover(expressionvariables(i)),F_expand,extendedvariables,monomtable,where,level+1,options,method,[],'concave');
140                    elseif isequal(properties.monotonicity,'exact')
141                        [F_expand,failure,cause] = expandrecursive(recover(expressionvariables(i)),F_expand,extendedvariables,monomtable,where,level+1,options,method,[],goal_vexity);
142                    else
143                        if options.allowmilp
144                            [F_expand,failure,cause] = expandrecursive(recover(expressionvariables(i)),F_expand,extendedvariables,monomtable,where,level+1,options,'milp',[],'milp');
145                        else
146                            failure = 1;
147                            cause = ['Monotonicity required at ' where ' at level ' num2str(level)];
148                        end
149                    end
150
151                end
152            end
153        end
154        if ~do_not_check_nonlinearity  & ~DUDE_ITS_A_GP & ~options.expandbilinear & ~options.allownonconvex
155            if isa(expression,'sdpvar')
156                if degree(expression)~=1 &~is(expression,'sigmonial')
157                    [Q,c,f,x,info] = quaddecomp(expression);
158                    if info
159                        failure = 1;
160                        cause = ['Polynomial expression  at ' where ' at level ' num2str(level)];
161                    else
162                        eigv = real(eig(Q));
163                        if ~all(diff(sign(eigv))==0)
164                            failure = 1;
165                            cause = ['Indefinite quadratic in ' where ' at level ' num2str(level)];
166                        else
167                            fail1 =  isequal(goal_vexity,'convex')  & all(eigv<=0) &  ~isequal(properties.monotonicity,'decreasing');
168                            fail2 =  isequal(goal_vexity,'convex')  & all(eigv>0)  &  ~isequal(properties.monotonicity,'increasing');
169                            fail3 =  isequal(goal_vexity,'concave') & all(eigv<=0) &  ~isequal(properties.monotonicity,'increasing');
170                            fail4 =  isequal(goal_vexity,'concave') & all(eigv>0)  &  ~isequal(properties.monotonicity,'decreasing');
171
172                            if fail1 | fail3
173                                failure = 1;
174                                cause = ['Concave quadratic encountered in ' where ' at level ' num2str(level)];
175                            elseif fail2 | fail4
176                                failure = 1;
177                                cause = ['Convex quadratic encountered in ' where ' at level ' num2str(level)];
178                            end
179                        end
180                    end
181                end
182            end
183        end
184        j = j+1;
185    end
186    F_expand = F_expand + F_graph;
187end
Note: See TracBrowser for help on using the repository browser.