www.gusucode.com > mbcdata 工具箱 matlab 源码程序 > mbcdata/@cgoptimrunner/private/pGradient.m

    function GFunc = pGradient(obj)
%PGRADIENT Calculate the gradient of expressions.
%
%   GFUNC = PGRADIENT(OBJ) generates a function handle to a function that
%   will calculates the gradients of the expressions in the provided
%   expression store.  The function handle must be called with the syntax G
%   = GFUNC(OBJ, X, Y, ITEMIND, ITEMDATA).  G is a cell array that contains
%   a cell array of expression gradients for each item.  X is a matrix of
%   free variable values, Y is a cell array of evaluation outputs, ITEMIND
%   is the items that expression gradients are required for and SDATA is
%   the cell array of expression value/input dependencies.

%   Copyright 2005-2009 The MathWorks, Inc.


% Create some indices for each free variable's values in the free value
% matrix.  This saves having to calculate it during every iteration.
FreeVarIndices = cell(1, length(obj.FreeVariableIndices));
End = 0;
for i = 1:length(FreeVarIndices)
    Start = End+1;
    End = End + obj.InputDataLengths(obj.FreeVariableIndices(i));
    FreeVarIndices{i} = Start:End;
end
GFunc = iCreateGradFunctions(FreeVarIndices);
end

function GFunc = iCreateGradFunctions(FreeVarIndices)

% Work vector storage for numjac: {Items, Fac, G}.  The first row is
% reserved for "All objectives" evaluation.  Additional lines are added as
% new combinations of ItemInd are asked for
WorkData = cell(1,3);
Threshold = [];
GFunc = @i_EvalGradient;

    function [G, FreeValIdx] = i_EvalGradient(obj, X, Y, ItemInd, ItemData, varargin)
        if isempty(Y)
            % No items to get gradient for, so just return an empty
            G = {};
            FreeValIdx = {};
            return
        end
        
        if nargin<4
            ItemInd = [];
        end
        ExprGroup = ItemData.ExpressionStore;
        
        % Find work data for these items
        idx = i_FindWorkData(ItemInd);
        
        % Generate basic sparsity pattern for the items requested.
        if isempty(ItemInd)
            S = vertcat(ItemData.ExpressionDeps{ItemData.ItemDataOffset+(1:length(ItemData.Items))});
        else
            S = vertcat(ItemData.ExpressionDeps{ItemInd});
        end
        
        [NRuns, NVars] = size(X);
        if size(S,1)==0
            % There are no expressions to get the gradient for
            Gmat = zeros([NRuns, 0, NVars]);
        else
            
            % Replicate sparsity pattern for each evaluation point we are
            % working at
            Sall = cell(1, NRuns);
            [Sall{:}] = deal(S);
            S = blkdiag(Sall{:});
            
            % Transform free variable data and evaluation outputs into a
            % suitable form
            Xnj = X.';
            Xnj = Xnj(:);
            Ynj = [Y{:}];
            Ynj = [Ynj{:}].';
            Ynj = Ynj(:);
            
            if isempty(Threshold)
                % just calculate threshold based on initial value
                Threshold = iCalcThreshold(obj,Xnj);
            end
            
            % Call numjac
            [BigG,WorkData{idx,2},WorkData{idx,3}] = ...
                numjac(@i_Evaluate, [], Xnj ,Ynj, ...
                Threshold, WorkData{idx,2}, 1, S, WorkData{idx,3}, ...
                obj, NRuns, NVars, ExprGroup, ItemData, ItemInd,varargin{:});
            
            
            % Transform the gradients into the appropriate groups to match the
            % original function evaluation format
            
            % Pull the blocks of gradients off the diagonal
            NFunVals = round(length(Ynj)./NRuns);
            Gmat = zeros([NRuns, NFunVals, NVars]);
            FunIdx = 1:NFunVals;
            VarIdx = 1:NVars;
            for n = 1:NRuns
                Gmat(n, :, :) = BigG(FunIdx, VarIdx);
                VarIdx = VarIdx + NVars;
                FunIdx = FunIdx + NFunVals;
            end
        end
        
        % Split this 3D Gmat across the second dimension.  Also pull out
        % only the relevant 3rd-dimension columns - ie only the free
        % variables that each objective item depends on
        G = cell(size(Y));
        SplitEnd = 0;
        if isempty(ItemInd)
            ItemInd = 1:length(Y);
        end
        FreeValIdx = cell(size(Y));
        for n = 1:length(Y)
            G{n} = cell(size(Y{n}));
            
            % Need to extract df/dy for variables in the order they are
            % reported by the item
            [IsFreeVar, FreeVarsIdx] = ismember(ItemData.InputIndices{ItemInd(n)}, obj.FreeVariableIndices);
            FreeValIdx{n} = [FreeVarIndices{FreeVarsIdx(IsFreeVar)}];
            for m = 1:length(Y{n})
                SplitStart = SplitEnd + 1;
                SplitEnd = SplitStart + size(Y{n}{m}, 2) - 1;
                G{n}{m} = Gmat(:, SplitStart:SplitEnd, FreeValIdx{n});
            end
        end
    end

    function idx = i_FindWorkData(ItemInd)
        if isempty(ItemInd)
            idx = 1;
        else
            idx = 0;
            for n = 2:size(WorkData,1)
                if isequal(WorkData{n,1}, ItemInd)
                    idx = n;
                    break
                end
            end
            if idx==0
                % Add a new entry
                WorkData = [WorkData; {ItemInd, [], []}];
                idx = size(WorkData,1);
            end
        end
    end

end

% Evaluation routine for the numjac call
function y = i_Evaluate(~, x, obj, NRuns, NVars, ExprGroup, ItemData, ItemInd, varargin)
NEvalPts = size(x, 2);
% Reconstruct free variable matrix from the vector x
x = reshape(x, NVars, NRuns*NEvalPts).';
Xeval = pCreateEvalInputs(obj, x);
if isempty(ItemInd)
    Y = pEvaluate(obj, ExprGroup, Xeval, ...
        ItemData.ItemDataOffset+(1:length(ItemData.Items)), varargin{:});
else
    Y = pEvaluate(obj, ExprGroup, Xeval, ItemInd, varargin{:});
end
y = [Y{:}];
y = [y{:}].';
NFunVals = size(y, 1);
y = reshape(y, NFunVals*NRuns, NEvalPts);
end


function Thresh = iCalcThreshold(obj,X)

Threshold = max(abs(double(X)*1e-8),1e-6);
TypicalY = double(X);
% set typical value to 1 where V is 0 otherwise numjac is poor
Index0 = abs(TypicalY)<1e-6;
if any(Index0)
    % use ranges
    [lb, ub] = getFreeVariableRanges(obj);
    TypicalY(Index0) = (lb(Index0)+ub(Index0))/2;
    Index0 = abs(TypicalY)<1e-6;
    if any(Index0)
        % use maximum bound
        MaxB =  max(abs(lb),abs(ub));
        TypicalY(Index0) = MaxB(Index0);
    end
end

Thresh = {Threshold,TypicalY};
end