[37] | 1 | function [F_expand,failure,cause] = expandrecursive(variable,F_expand,extendedvariables,monomtable,where,level,options,method,extstruct,goal_vexity) |
---|
| 2 | global DUDE_ITS_A_GP ALREADY_MODELLED REMOVE_THESE_IN_THE_END MARKER_VARIABLES OPERATOR_IN_POLYNOM |
---|
| 3 | cause = ''; |
---|
| 4 | failure = 0; |
---|
| 5 | |
---|
| 6 | if ~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; |
---|
| 187 | end |
---|