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

    function script = getJacobScript(m)
%X2FX  regression X matrix for xregcubic
%
% FX= x2fx(m,X)

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

N = order(m);
maxInteract = get(m, 'maxInteract');
numOfCoeffs = size(double(m), 1);
reorderM = reorder(m);

N(maxInteract+1:end)=0;

% recusive algorithm for Jacobian
FX = repmat({''}, 1, numOfCoeffs);
FX{1} = '1';
FX = recursiveJacob(FX,1,1, N, maxInteract, FX{1}, 1, reorderM);

% add higher order terms
FX = iAddHigherOrderTerms(m, FX, reorderM);

% add the begining part "FX(i) = ", to the FX values
for j=1:size(FX,2)
    FX{j} = sprintf('FX(%d) = %s;', j,FX{j});
end
FX = sprintf('%s\n', FX{:});

% combine header with help
header = {
    'function FX = fcn(X)'
    '%#codegen'
    '% Jacobian calculation for a polynomial model'
    ''
    sprintf('%% Autogenerated by the Model-Based Calibration Toolbox on %s.', datestr(now));
    ''
    sprintf('FX = zeros(%d,1);', numOfCoeffs);
    };
header = sprintf('%s\n', header{:});
% add code to calculate the Jacobian
script = sprintf('%s%s', header, FX);


function FX = iAddHigherOrderTerms(m, FX, reorderM)
N = order(m);
maxInteract = get(m, 'maxInteract');

% non interactive terms
nx = sum(N(maxInteract+1:end));
if nx>0
    p = size(FX,2)-nx;
    % Xm = {'X(1)', 'X(2)', ..., 'X(n)'}, text used to the ith element of X 
    Xm = arrayfun(@(i)sprintf('X(%d)',reorderM(i)), 1:N(maxInteract+1), 'UniformOutput', 0);
    % Xi is the intermediate value
    % initial value: Xi(i) = Xm(i)^maxInteract
    Xi = cellfun(@(x)sprintf('%s^%d',x,maxInteract), Xm, 'UniformOutput', 0);
    for i = maxInteract+1:length(N)
        % Xi = Xi.*Xm;
        Xi = cellfun(@(x,y)[x,'*',y], Xm, Xi, 'UniformOutput', 0);
        
        for j=1:N(i)
            p = p+1;
            FX(p) = Xi(j);
            Xi{j} = sprintf('FX(%d)',p);
        end
    end
end


function [FX, currpos] = recursiveJacob(FX,lvl,start, N, maxInteract, prefix, currpos, reorderM)

for i=start:N(lvl)
    currpos = currpos+1;
    if strcmp(prefix, '1') % can skip multiplication when mulitplying by 1
        FX{currpos} = sprintf('X(%d)', reorderM(i));
    else
        FX{currpos} = sprintf('%s*X(%d)', prefix, reorderM(i));
    end
    newPrefix = sprintf('FX(%d)', currpos);
    if lvl < maxInteract
        [FX, currpos] = recursiveJacob(FX,lvl+1,i, N, maxInteract, newPrefix, currpos, reorderM);
    end
end