expandrecursive.m

matlab波形优化算法经常要用到的matlab toolbox工具箱:yalmip

function [F_expand,failure,cause] = expandrecursive(variable,F_expand,extendedvariables,monomtable,variabletype,where,level,options,method,extstruct,goal_vexity,allExtStruct)

% Some crazy stuff to help out when hand;ling nonocnvex and other
% non-standard problems
global DUDE_ITS_A_GP ALREADY_MODELLED ALREADY_MODELLED_INDEX REMOVE_THESE_IN_THE_END MARKER_VARIABLES OPERATOR_IN_POLYNOM

% Assume no failure due to failed model or convexity propagation
cause = '';
failure = 0;

% We only expand this variable if it hasn't already been modelled in
% another constraint, and that model is sufficicent for the current case.
if  ~alreadydone(getvariables(variable),method,goal_vexity)

    % Request model for this variable
    ext_index = find(getvariables(variable) == extendedvariables);

    % [F_graph,properties,arguments,fcn] = model(variable,method,options,extstruct);
    [properties,F_graph,arguments,fcn] = model(variable,method,options,allExtStruct(ext_index));

    % Bit of a hack (CPLEX sucks on some trivial problems, so we have to track
    % these simple models and fix them in the end). Basically, when a
    % constraint such as ismember(x,B) is used, YALMIP introduces a "marker
    % variable" and returns the constraint marker == 1. This constraint and
    % variable should be removed in the end, it is only used to ensure that
    % the operator code is run. CPLEX can fail on this idiot constraints. 
    if ~isempty(properties)
        if isfield(properties,'extra')
            if strcmp(properties.extra,'marker')
                MARKER_VARIABLES = [MARKER_VARIABLES getvariables(variable)];
            end
        end
    end

    % If we are running a graph representation, we have to make sure that
    % we really are satisfying the convexity constraints    
    if isequal(method,'graph')

        failure = ~isequal(properties.convexity,goal_vexity);

        if failure & ~options.allownonconvex
            failure = 1;
            cause = ['Expected ' goal_vexity ' function in '  where ' at level ' num2str(level)];
            return
        end
       
        if failure & any(strcmpi(properties.model,{'callback','exact','integer'}))
            % This operator is guaranteed to return the associated value
            % (i.e. it is not modelled by graphs) and can thus be used in
            % the model. However, we have to give up convexity propagation
            % since it didn't satisfy the convexity requirements
            failure = 0;
            method = 'exact';
        else
            if failure & options.allowmilp
                [properties,F_graph,arguments] = model(variable,'exact',options);              
                if isempty(F_graph)
                    cause = [cause ', exact model not available.'];
                    return
                else
                    failure = 0;                   
                    method = properties.model;
                end
            elseif failure
                cause = ['Expected ' 'goal_vexity' ' function in '  where ' at level ' num2str(level)];
            end
        end
    else
        if isempty(F_graph)
            cause = ['Model not available in ' where ' at level ' num2str(level)];
            failure = 1;
            return
        else
            failure = 0;
            method = 'exact';
        end
    end

    % We save variable models for future use
   if failure == 0
   done = save_model_expansion(method,goal_vexity,F_graph,properties);
   if done
       return
   end
   end

    % Now we might have to recurse
    if max(size(arguments))>1
        arguments = reshape(arguments,prod(size(arguments)),1);
    end

    [ix,jx,kx] = find(monomtable(getvariables(variable),:));
    if ~isempty(jx) % Bug in 6.1
        if any(kx>1)
            OPERATOR_IN_POLYNOM = [OPERATOR_IN_POLYNOM extendedvariables(jx(find(kx>1)))];
        end
    end

    % Small pre-processing to speed-up large-scale problems (subsref sloooow)
    % with only linear arguments (such as norm(Ax-b) problems)
    if isa(arguments,'sdpvar')
        argvars = getvariables(arguments);
        if max(argvars)>length(variabletype)
            variabletype = yalmip('variabletype');
        end
        do_not_check_nonlinearity = ~any(variabletype(argvars));
        if do_not_check_nonlinearity
            allvariables = argvars;
            fullbasis = getbase(arguments);
            fullbasis = fullbasis(:,2:end);
            fullbasis_transpose = fullbasis';
        end
    else
        do_not_check_nonlinearity = 0;
    end

    j = 1;
    % Ok, here goes the actual recursive code
    while j<=length(arguments) & ~failure

        if do_not_check_nonlinearity            
            usedvariables = find(fullbasis_transpose(:,j));
            expressionvariables = allvariables(usedvariables);
        else
            expression = arguments(j);
            expressionvariables = unique([depends(expression) getvariables(expression)]);
        end
        index_in_expression = find(ismembc(expressionvariables,extendedvariables));

        if ~isempty(index_in_expression)
            for i = index_in_expression
                if ~milpalreadydone(expressionvariables(i))
                    if do_not_check_nonlinearity                       
                        basis = fullbasis(j,usedvariables((i)));
                    else
                        basis = getbasematrix(expression,expressionvariables(i));
                    end
                    go_convex1 = (basis > 0) &  isequal(goal_vexity,'convex') & isequal(properties.monotonicity,'increasing');
                    go_convex2 = (basis <= 0) &  isequal(goal_vexity,'convex') & isequal(properties.monotonicity,'decreasing');
                    go_convex3 = (basis <= 0) &  isequal(goal_vexity,'concave') & isequal(properties.monotonicity,'increasing');
                    go_convex4 = (basis > 0) &  isequal(goal_vexity,'concave') & isequal(properties.monotonicity,'decreasing');

                    go_concave1 = (basis > 0) &  isequal(goal_vexity,'convex') & isequal(properties.monotonicity,'decreasing');
                    go_concave2 = (basis <= 0) &  isequal(goal_vexity,'convex') & isequal(properties.monotonicity,'increasing');
                    go_concave3 = (basis <= 0) &  isequal(goal_vexity,'concave') & isequal(properties.monotonicity,'decreasing');
                    go_concave4 = (basis > 0) &  isequal(goal_vexity,'concave') & isequal(properties.monotonicity,'increasing');

                    if go_convex1 | go_convex2 | go_convex3 | go_convex4
                        [F_expand,failure,cause] = expandrecursive(recover(expressionvariables(i)),F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,method,[],'convex',allExtStruct);
                    elseif go_concave1 | go_concave2 | go_concave3 | go_concave4
                        [F_expand,failure,cause] = expandrecursive(recover(expressionvariables(i)),F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,method,[],'concave',allExtStruct);
                    elseif isequal(properties.monotonicity,'exact')
                        [F_expand,failure,cause] = expandrecursive(recover(expressionvariables(i)),F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,method,[],goal_vexity,allExtStruct);
                    else
                        if options.allownonconvex%milp
                            [F_expand,failure,cause] = expandrecursive(recover(expressionvariables(i)),F_expand,extendedvariables,monomtable,variabletype,where,level+1,options,'exact',[],'none',allExtStruct);
                        else
                            failure = 1;
                            cause = ['Monotonicity required at ' where ' at level ' num2str(level)];
                        end
                    end

                end
            end
        end
        if ~do_not_check_nonlinearity  & ~DUDE_ITS_A_GP & ~options.expandbilinear & ~options.allownonconvex
            if isa(expression,'sdpvar')
                if degree(expression)~=1 &~is(expression,'sigmonial')
                    [Q,c,f,x,info] = quaddecomp(expression);
                    if info
                        failure = 1;
                        cause = ['Polynomial expression  at ' where ' at level ' num2str(level)];
                    else
                        eigv = real(eig(Q));
                        if ~all(diff(sign(eigv))==0)
                            failure = 1;
                            cause = ['Indefinite quadratic in ' where ' at level ' num2str(level)];
                        else
                            fail1 =  isequal(goal_vexity,'convex')  & all(eigv<=0) &  ~isequal(properties.monotonicity,'decreasing');
                            fail2 =  isequal(goal_vexity,'convex')  & all(eigv>0)  &  ~isequal(properties.monotonicity,'increasing');
                            fail3 =  isequal(goal_vexity,'concave') & all(eigv<=0) &  ~isequal(properties.monotonicity,'increasing');
                            fail4 =  isequal(goal_vexity,'concave') & all(eigv>0)  &  ~isequal(properties.monotonicity,'decreasing');

                            if fail1 | fail3
                                failure = 1;
                                cause = ['Concave quadratic encountered in ' where ' at level ' num2str(level)];
                            elseif fail2 | fail4
                                failure = 1;
                                cause = ['Convex quadratic encountered in ' where ' at level ' num2str(level)];
                            end
                        end
                    end
                end
            end
        end
        j = j+1;
    end
    F_expand = F_expand + F_graph;
end