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