optimShape

OPTIMSHAPE Finds the optimal shape parameter of a RBF. Also returns the optimal
interpolant
 
 Given known input data, it uses a percentage of this data 
 (1-validationPercent) to create an interpolant. Then, it uses the rest of
 the data to validate the interpolant accuracy (using Sum of Squared
 Differences). This function is used in a non-linear optimization to seek
 for the optimal 'e' (shape parameter of some RBF)

0001 function [e, res, testInd, validationInd] = optimShape(x, y, z, rbfType, validationPercent, eInit, options)
0002 %OPTIMSHAPE Finds the optimal shape parameter of a RBF. Also returns the optimal
0003 %interpolant
0004 %
0005 % Given known input data, it uses a percentage of this data
0006 % (1-validationPercent) to create an interpolant. Then, it uses the rest of
0007 % the data to validate the interpolant accuracy (using Sum of Squared
0008 % Differences). This function is used in a non-linear optimization to seek
0009 % for the optimal 'e' (shape parameter of some RBF)
0010 %
0011 
0012 if strcmpi(rbfType, 'linear') || ...
0013     strcmpi(rbfType, 'cubic') || ...
0014     strcmpi(rbfType, 'quintic') || ...
0015     strcmpi(rbfType, 'thinplate') || ...
0016     strcmpi(rbfType, 'green') 
0017     error('%s RBF does not require a shape parameter!', rbfType);
0018 end
0019 
0020 if numel(x) ~= numel(y) || numel(x) ~= numel(z)
0021     error('The number of elements in x, y, z must be the same!');
0022 end
0023 
0024 if validationPercent <= 0 || validationPercent >= 1
0025     errro('validationPercent must be specified as a number between 0 and 1 (non-inclusive)');
0026 end
0027 
0028 % Separate the reference data into test and validation
0029 numValid = ceil(numel(x)*validationPercent);
0030 rp = randperm(length(x));
0031 validationInd = rp(1:numValid);
0032 testInd = rp(numValid+1:end);
0033 
0034 % Test data
0035 xt = x(testInd);
0036 yt = y(testInd);
0037 zt = z(testInd);
0038 
0039 % Validation data
0040 xv = x(validationInd);
0041 yv = y(validationInd);
0042 zv = z(validationInd);
0043 
0044 % Optimize
0045 objectiveFun = @(e) residuals(xt, yt, zt, xv, yv, zv, options, rbfType, e);
0046           
0047 % We turn off the MATLAB:nearlySingularMatrix warning, as it will happen a
0048 % lot during optimization
0049 warning('off', 'MATLAB:nearlySingularMatrix');
0050 
0051 % Minimize
0052 % optimOptions = optimset('MaxFunEvals', 10000000, 'TolFun', 1e-10, 'TolX', 1e-10);
0053 optimOptions = optimset('MaxFunEvals', 10000000);
0054 [e, res, exitflag, output] = fminsearch(objectiveFun, eInit, optimOptions);
0055 
0056 % Turn on the warning again
0057 warning('on', 'MATLAB:nearlySingularMatrix');
0058 end
0059 
0060 
0061 function res = residuals(xt, yt, zt, xv, yv, zv, options, rbfType, e)
0062 % Input:
0063 %   - xt, yt, zt: test values
0064 %   - xv, yv, zv: validation values
0065 %
0066 % Output:
0067 %   - res = summed residuals
0068 
0069 rbfInterp = RBFInterpolant(xt, yt, zt, 'DistanceType', options.DistanceType, ...
0070                                        'PolynomialDegree', options.RBF.(rbfType).PolynomialDegree, ...
0071                                        'RBF', rbfType, ...
0072                                        'RBFEpsilon', e, ... 
0073                                        'Smooth', options.RBF.(rbfType).Smooth, ...
0074                                        'Regularization', options.RBF.(rbfType).Regularization);
0075 zi = rbfInterp.interpolate(xv, yv);
0076 
0077 res = sum(abs(zv-zi));
0078 
0079 end