www.gusucode.com > mbcmodels 工具箱 matlab 源码程序 > mbcmodels/@xreginterprbf/interpolate.m

    function [om, OK] = interpolate( m )
%XREGINTERPRBF/INTERPOLATE   Find rbf coefficients by interpolation
%  [OM,OK] = INTERPOLATE(M)

%  Copyright 2000-2009 The MathWorks, Inc. and Ford Global Technologies, Inc.


om = contextimplementation( xregoptmgr, m, @i_interpolate, [], ...
    'Interpolation', @interpolate );

% fit parameters
om = AddOption( om, 'PolyFromKernel', ...
    1, 'boolean', ...
    'Polynomial from kernel', false );

om = AddOption( om, 'CoincidentStrategy', ...
    'Mean', 'Maximum|Minimum|Mean', ...
    'Coincident Node Strategy', true );

om = AddOption( om, 'Algorithm', ...
    'Direct', 'Direct|GMRES|BICG|CGS|QMR', ...
    'Algorithm', true );

om = AddOption( om, 'Tolerance', ...
    1e-6, {'numeric',[0, Inf]}, ...
    'Tolerance', true );

om = AddOption( om, 'MaxIt', ...
    20, {'int', [1, Inf]}, ...
    'Maximum number of iterations', true );

om = AddOption( om, 'cost', ...
    0, {'numeric', [-Inf, Inf]}, ...
    [], false );

OK = 1;

return

% -------------------------------------------------------------------------|
function [m, cost,OK ] = i_interpolate( m, om, x0, x, y, varargin )
%  Inputs:		m   xreinterprbf object
%               om  xregoptmgr
%               x0  starting values (not used)
%				x   matrix of data points 
%				y   target values
%
% Outputs:     m    new rbf object
%              cost log10GCV

% Get the Coincident Node Strategy from the opt-mgr
% -- Note that the code historically had two speelings of this property.
%    Hence we try the one we state above where we set up the opt-mgr but
%    just in case that doesn't work we try the alternative spelling.
try
    coincidentStrategy = get( om, 'CoincidentStrategy' );
catch e
    if strcmpi( e.identifier, 'mbc:xregoptmgr:InvalidPropertyName' ),
        % this hides a typo
        OldPropName = char([67,111,110,105,110,99,105,100,101,110,116, 83,116,114 97,116,101,103,121]);
        coincidentStrategy = get( om, OldPropName );
    else
        rethrow( e );
    end
end

% should the polynomial degree be set by the kernel
if get( om, 'PolyFromKernel' ),
    m = check_degree( m, 'Hard' );
end

% if using a thin-plate, linear or cubic kernel, set the width to 1.0
if any( strcmpi( get( m, 'kernel' ), {'thinplate', 'linearrbf', 'cubicrbf'} ) ),
    m = setrbfpart( m, 'width', 1.0 );
end

% set the centers of the rbf to the input nodes
[centers, ~, J] = unique( x, 'rows' ); % remove identical nodes
ncenters = size( centers, 1 );         % number of rbf centers
if ncenters == size( x, 1 ),
    % all nodes are centers,
    Y = y;
    centers = x;
else
    switch lower( coincidentStrategy ),
    case 'maximum'
        alg = @max;
    case 'minimum'
        alg = @min;
    case 'mean'
        alg = @mean;
    otherwise
        error(message('mbc:xreginterprbf:InvalidState', coincidentStrategy));
    end
    Y = zeros( ncenters, 1 );
    for i = 1:ncenters,
        Y(i) = feval( alg, y(J==i) );
    end
end

poly = get( m, 'linearmodpart' );
poly_dim = size( poly, 1 );

% make sure the status of the rbf centers is ok.
current_status = get( m, 'status' );
m = setrbfpart( m, 'centers', centers );
m = setstatus( m, poly_dim+(1:ncenters), 1 ); % 1 == always
if poly_dim+ncenters < size( current_status, 1 ),
    % the new center list is shorter than the old.
    ind = (poly_dim+ncenters+1):size(current_status,1);
    m = setstatus( m, ind, 2 ); % 2 == never
end

% make sure the status of the linearmodpart is correct
nt=size(get( m, 'linearmodpart' ),1)+ncenters;
Tin = true( nt, 1 );
Tin(1:poly_dim)= Terms(poly);

% setup the interpolation system
FX = CalcJacob( m, centers );
n = size(FX,2) - ncenters; % dimension of polynomial space
if n > 0,
    P = FX(:,1:n);
    A = [ FX; zeros( n ),  P' ];
    b = [Y; zeros(n,1)];
else
    A = FX;
    b = Y;
end

% solve the interpolation system
maxit = get( om, 'MaxIt' );
tol = get( om, 'Tolerance' );
switch get( om, 'Algorithm' ),
case 'Direct',
    if isempty(A)
        flag = true;
        coeffs = zeros(sum(Tin),1);
    else
        [coeffs,r] = linsolve(A,b);
        % use inverse condition number to determine rank deficiency
        flag = r < length(A)*eps( max(abs(A(:))) );
    end
case 'GMRES',
    % maximum iterations <= size(A,1)
    maxit= min(maxit,size(A,1));
    [coeffs, flag] = gmres(A,b,[],tol,maxit);
case 'BICG',
    [coeffs, flag] = bicg(A,b,tol,maxit);
case 'CGS',
    [coeffs, flag] = cgs(A,b,tol,maxit);
case 'QMR',
    [coeffs, flag] = qmr(A,b,tol,maxit);
otherwise
    error(message('mbc:xreginterprbf:InvalidState1', get( om, 'Algorithm' )));
end

beta= zeros(nt,1);
beta(Tin)= coeffs;

% set model coefficients
m = update( m, beta );

% setup the cost statistic
cost = 0;
setFitOpt( m, 'cost', cost ); 
OK = ~flag;
%--------------------------------------------------------------------------
% EOF
%--------------------------------------------------------------------------