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

    function varargout = pRBFExtrapSurfaceGen(axisbp, pts, x, action, SurfaceGenData)
%PRBFEXTRAPSURFACEGEN RBF Surface Generation function
%
%   PRBFEXTRAPSURFACEGEN uses RBF extrapolation to generate a surface
%   utilising the extrapolation routines in CGMATHSOBJECT. This function is
%   called by CGTABGRADCONSRAINT using the syntaxes below.
%
%   SURFACEGENDATA = PRBFEXTRAPSURFACEGEN(AXISBP, PTS, X, 'FIXEDDATA')
%   calculates the interpolation matrices and the duplicate transform.
%   These matrices only depend on the fixed variables and will be stored in
%   the working store of the table gradient constraint. 
%
%   [Y, YGRAD] = PRBFEXTRAPSURFACEGEN(AXISBP, PTS, X, 'EVALUATE',
%   SURFACEGENDATA) returns an extrapolation of X from PTS to the axis
%   breakpoints AXISBP. The gradient of this function is returned in YGRAD.
%   Note that the axis variables are scaled onto [0, 1]. 
%
%   [Y, YGRAD] = PRBFEXTRAPSURFACEGEN(AXISBP, PTS, X, 'EVALUATE',
%   SURFACEGENDATA) returns a extrapolation of multiple values of X in the
%   case where X is a NPTS-by-NRUN matrix. In this case, the axis values,
%   PTS, are the same for each column of X. Y is a NTAB-by-NRUN matrix of
%   extrapolations of X and YGRAD is a NTAB-by-NPTS-by-NRUN matrix of
%   gradient evaluations.
%
%   See also PRBFSURFACEGEN

%   Copyright 2006-2015 The MathWorks, Inc.

switch lower(action)
    case 'fixeddata'
        varargout{1} = i_FixedData(axisbp, pts);
    case 'evaluate'
        [varargout{1:nargout}] = i_Evaluate(axisbp, pts, x, SurfaceGenData);
end

end
%--------------------------------------------------------------------------
function SurfaceGenData = i_FixedData(~, ~)
%--------------------------------------------------------------------------

% No need to store any data for this function
SurfaceGenData = [];

end
%--------------------------------------------------------------------------
function varargout = i_Evaluate(axisbp, pts, x, ~)
%--------------------------------------------------------------------------

% Initialise return
NTABPTS = prod(cellfun('length', axisbp));
NRUN = size(x, 2);
varargout{1} = zeros(NTABPTS, NRUN);
if nargout > 1
    varargout{2} = zeros([NTABPTS, size(x)]);
end

% Loop over specified runs
for i = 1:size(x, 2)

    % Call extrapolate fill routine
    NAXES = length(axisbp);
    varargout{1}(:, i) = n_FillTable(x(:, i));

    % Calculate the gradient of the extrapolation
    if nargout > 1
        varargout{2}(:, :, i) = n_CalcGradient(x(:, i));
    end
    
end

    %----------------------------------------------------------------------
    function newtabdata = n_FillTable(x)
    %----------------------------------------------------------------------

        % Generate fill data for extrapolation routine
        filldata = [pts, x];

        % Temporary copy from filling routine in CAGE Optimization. If this
        % test is successful, I'll look into commonizing this routine.

        switch NAXES
            case 2
                % 2D table

                % Swap axes order for old filling routine
                axes = axisbp([2 1]);
                filldata = filldata(:, [2 1 3]);

                % Do the extrapolation
                newtabdata = eval(cgmathsobject,'extrapolate_RBF', ...
                    filldata(:,1), filldata(:, 2), filldata(:,3), axes{1}, axes{2}); %#ok<GTARG>

            case 1
                % 1D Table

                % Get fill data
                [xpts, xi] = unique(filldata(:,1)');
                ypts = filldata(xi, 2);

                % Do the extrapolation
                newtabdata = eval(cgmathsobject, 'extinterp1', xpts, ...
                    ypts, axisbp{1}(:)); %#ok<GTARG>

        end
        
        % Return a column vector
        newtabdata = newtabdata(:);

    end

    %----------------------------------------------------------------------
    function ygrad = n_CalcGradient(x)
    %----------------------------------------------------------------------

        % Use the finite differences code from Optimization Toolbox
         
        % Inputs for finite differences
        objfun = @n_FillTable;
        lb = [];
        ub = [];
        fCurrent = n_FillTable(x);
        numberOfVariables = numel(x); 
        variables = 1:numberOfVariables;
        gradf = zeros(length(fCurrent), numberOfVariables);
        finDiffFlags.chkFunEval = false;
        finDiffFlags.fwdFinDiff = true;
        finDiffFlags.isGrad = false;
        finDiffFlags.hasLBs = false(numberOfVariables,1);
        finDiffFlags.hasUBs = false(numberOfVariables,1);
        finDiffFlags.scaleObjConstr = false;
        finDiffOpts.FinDiffType = 'forward';
        finDiffOpts.FinDiffRelStep = sqrt(eps)*ones(numberOfVariables,1);
        finDiffOpts.TypicalX = ones(size(x(:)));
        finDiffOpts.DiffMinChange = 1e-4;
        finDiffOpts.DiffMaxChange = 1e-1;
        % sizes structure
        sizes.xShape = size(x);
        sizes.nVar = numberOfVariables;
        sizes.mNonlinIneq = 0;
        sizes.mNonlinEq = 0;
        
        % Call finite differences
        ygrad = finitedifferences(x(:), objfun, [], lb, ub, fCurrent, [], ...
            [], variables, finDiffOpts, sizes, gradf, [], [], finDiffFlags, []);
        
    end
end