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 %--------------------------------------------------------------------------