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

    function script = getScript(m)
%getScript Creates the script to evaluate an RBF model
% Generates a script for rbf evaluation and Jacobian.
%
% s = getScript(M) returns a script as a char array that
% evaluates an RBF model M.

%  Copyright 2013 The MathWorks, Inc. and Ford Global Technologies, Inc.

kernelName = getkernelstring(m);
width = m.width;

% 2-D widths case = [#centers x #factors]
twoDCase = size(width, 1) > 1 && size(width,2)==nfactors(m);
% isWidthPerDim = width is a vector or the 2-D Case
isWidthPerDim = (numel(width)==nfactors(m)) || twoDCase;
isemptyW = isempty(width);

headerText = generateHeader(twoDCase, kernelName);

kernelText = KernelText(kernelName, isWidthPerDim, isemptyW);

sqDistTextCell = RBFDistanceText(isWidthPerDim, isemptyW, kernelName);
sqDistText = sprintf('%s\n', sqDistTextCell{:});

script = sprintf('%s\n\n%s\n\n%s', headerText, kernelText, sqDistText);


function headerText = generateHeader(twoDCase, kernelName)

fcnTextCell = {
    'function J = fcn(X, center, width)'
    '%#codegen'
    ''
    sprintf('% RBF evaluation for %s kernel', kernelName)
    '% X      = evaluation points [inputs x 1]'
    '% center = RBF centers [inputs x ypoints]'
    '% width  = width as a value (1x1), vector (1 x ypoints) or matrix (inputs x ypoints)'
    '% Output J = jacobian, can be combined with Beta to evaluate the RBF'
    ''
    sprintf('%% Autogenerated by the Model-Based Calibration Toolbox on %s.', datestr(now));
    'nC = size(center, 1);'
    'd = size(center, 2);'
    ''
    '%  Set up output, jacobian case'
    'J = zeros(nC,1);'
    ''
    '% Main Computation loop'
    'for i=1:nC'
    callToKernel(twoDCase)
    'end'};
headerText = sprintf('%s\n', fcnTextCell{:});


function text = callToKernel(twoDCase)
if twoDCase
    text = '    J(i) = rbfKernel(d, X, center(i,:), width(i,:));';
else
    text = '    J(i) = rbfKernel(d, X, center(i,:), width);';
end


function text = KernelText(kernelName, isWidthPerDim, isemptyW)
switch kernelName
    case 'distance'
        o = RBFDistanceText(isWidthPerDim, isemptyW, kernelName);
    case 'gaussian'
		o = RBFGaussianText();
    case 'multiquadric'
        o = RBFMultiquadricText(isWidthPerDim, isemptyW);
    case 'recmultiquadric'
		o = RBFRecMultiquadricText(isWidthPerDim, isemptyW);
    case 'linear'
		o = RBFLinearText();
    case 'linearrbf'
		o = RBFLinearText();
    case 'thinplate'
		o = RBFThinPlateText();
    case 'cubic'
		o = RBFCubicText();
    case 'cubicrbf'
		o = RBFCubicText();
    case 'logistic'
        o = RBFLogisticText();
    case 'logisticrbf'
		o = RBFLogisticText();
    case 'wendland0'
        o = RBFWendland0Text();
    case 'wendland2'
		o = RBFWendland2Text();
    case 'wendland4'
        o = RBFWendland4Text();
    case 'wendland6'
        o = RBFWendland6Text();
    otherwise
        error(message('mbc:xregrbf:InvalidModel'));
end
text = sprintf('%s\n', o{:});


function text = RBFDistanceText(isWidthPerDim, isemptyW, kernelName)

if any(strcmp(kernelName, {'multiquadric', 'recmultiquadric'})) ...
        && ~isWidthPerDim
    % special case for multiquadratic kernel the weight might be set to []
    % when the rbfDistanceSq function is called
    isemptyW = true;
end
    
text = {
    'function d2 = rbfDistanceSq(n, x, y, w)'
    '% Distance Squared'
    'd2 = 0;'};

if isWidthPerDim
    textS = {
        'for i=1:n '
        '    d2 = d2 + (x(i) - y(i))^2/w(i);'
        'end'};
    text = [text;textS];
else
    textS = {
        'for i=1:n'
        '    d2 = d2 + (x(i) - y(i))^2;'
        'end'};
    text = [text;textS];
    if ~isemptyW
        textS = {'d2 = d2/w;'};
        text = [text;textS];
    end
end


function text = RBFGaussianText()
text = {
    'function o = rbfKernel(n, x, y, w)'
    '% Gaussian'
    '% phi( r ) = exp( -r^2 )'
    'o = exp(-rbfDistanceSq(n, x, y, w));'};


function text = RBFMultiquadricText(isWidthPerDim, isemptyW, useActualName)
if nargin < 3 || ~useActualName
    fName = 'rbfKernel';
else
    fName = 'rbfMultiquadric';
end
text = {
    sprintf('function o = %s(n, x, y, w)', fName)
    '% Multiquadric'
    '%  phi(r) = sqrt((r/w)^2+1) if we have a width per dimension'
    '%  phi(r) = sqrt(r^2+w^2)   otherwise'};

if isWidthPerDim
    text(end+1) = {'o = sqrt(rbfDistanceSq(n, x, y, w)+1);'};
elseif ~isemptyW
    text(end+1) = {'o = sqrt(rbfDistanceSq(n, x, y, [])+w);'};
else
    text(end+1) = {'o = sqrt(rbfDistanceSq(n, x, y, [])+1);'};
end


function text = RBFRecMultiquadricText(isWidthPerDim, isemptyW)
textRecMultiquadratic = {
    'function o = rbfKernel(n, x, y, w)'
    '% Reciprocal Multiquadric'
    '%  phi(r) = 1/Multiquadric(r)'
    'o = 1/rbfMultiquadric(n, x, y, w);'};

textMultiquadratic = RBFMultiquadricText(isWidthPerDim, isemptyW, true);

text = [textRecMultiquadratic;{'';''};textMultiquadratic];


function text = RBFLinearText()
text = {
    'function o = rbfKernel(n, x, y, w)'
    '% Linear'
    '%  phi( r ) = r'
    'o = sqrt(rbfDistanceSq(n, x, y, w));'};


function text = RBFThinPlateText()
text = {
    'function o = rbfKernel(n, x, y, w)'
    '% Thin-plate spline'
    '%  phi(r) = r^2 log(r)'
    'r2 = rbfDistanceSq(n, x, y, w);'
    'if r2>0'
    '    o = 0.5*r2*log(r2);'
    'else'
    '    % log(0) case'
    '    o = 0;'
    'end'};


function text = RBFCubicText()
text = {
    'function o = rbfKernel(n, x, y, w)'
    '% Cubic'
    '%  phi(r) = r^3'
    'r = sqrt(rbfDistanceSq(n, x, y, w));'
    'o = r*r*r;'};


function text = RBFLogisticText()
text = {
    'function o = rbfKernel(n, x, y, w)'
    '% Logistic'
    '%  phi(r) = 1/(1 + exp(r/w))'
    'o = 1/(1+exp(sqrt(rbfDistanceSq(n, x, y, w))));'};


function text = RBFWendland0Text()
text = {
    'function phi = rbfKernel(n, x, y, w)'
    '% Wendland, continuity = 0'
    '%  phi(r) = (1-r)^(n/2+1)     0.0 <= r <= 1.0'
    '%  0                            1.0 <  r.'
    'r2 = rbfDistanceSq(n, x, y, w);'
    'if (r2>1)'
    '    phi = 0;'
    'else'
    '    pow = floor(n/2)+1;'
    '    phi = (1-sqrt(r2))^pow; '
    'end'};


function text = RBFWendland2Text()
text = {
    'function phi = rbfKernel(n, x, y, w)'
    '% Wendland, continuity = 2'
    '%  ell = n/2 + 2 '
    '%  phi(r) = ((ell+1)*r+1)*(1-r)^(ell+1),     0.0 <= r <= 1.0'
    '%           0,                               1.0 <  r.'
    'ell = floor(n/2)+2;'
    'r   = rbfDistanceSq(n, x, y, w);'
    'if r>1'
    '    phi = 0;'
    'else'
    '    r = sqrt(r);'
    '    phi = (ell+1)*r+1;'
    '    omr = (1-r)^(ell+1); % o(ne) m(inus) r'
    '    phi = phi*omr;'
    'end'};


function text = RBFWendland4Text()
text = {
    'function phi = rbfKernel(n, x, y, w)'
    '% Wendland, continuity = 4'
    '%  ell = n/2 + 3 '
    '%  phi( r ) = ((ell^2 + 4*ell + 3)*r^2 + (3*ell+6) * r + 3) * (1 - r)^(ell+2),     '
    '%                         0.0 <= r <= 1.0'
    '%             0,          1.0 <  r.'
    'ell = floor(n/2)+3;'
    'r = rbfDistanceSq(n, x, y, w);'
    'if r>1'
    '    phi = 0;'
    'else'
    '    r = sqrt(r);'
    '    phi = ((ell*(ell+4)+3)*r+(3*ell+6))*r+3;'
    '    omr = (1-r)^(ell+2); %  o(ne) m(inus) r'
    '    phi = phi*omr;'
    'end'};


function text = RBFWendland6Text()
text = {
    'function phi = rbfKernel(n, x, y, w)'
    '% Wendland, continuity = 6'
    '%  ell = n/2 + 4'
    '%  phi( r ) = 0,           1.0 <  r.'
    '%           = [(ell^3 + 9*ell^2 + 23*ell + 15)r^3 + (6*ell^2 + 36*ell + 45)*r^2 '
    '%                + (15*ell+45)r + 15] * (1 - r)^(ell+3),     0.0 <= r <= 1.0'
    'ell = floor(n/2)+4;'
    'r = rbfDistanceSq(n, x, y, w);'
    'if r>1'
    '    phi = 0;'
    'else'
    '    r = sqrt(r);'
    '    phi = (((((ell+9)*ell+23)*ell+15)*r+(6*ell+36)*ell+45)*r+15*ell+45)*r+15;'
    '    omr = (1-r)^(ell+3); % o(ne) m(inus) r'
    '    phi = phi*omr;'
    'end'};