function [S, R] = sensitivity(A,Aparamkey);% sensitivity.m Use finite-differences to compute sensitivity%               of a naemodel object with respect to a parameter.%%               Called as%%               [S, R] = sensitivity(A,Aparamkey);% Copyright (c) by Raymond A. Adomaitis, 1998-2005% Based on the finite-difference algorithm for computing% jacobian arrays (jacobian.m)% check if residual in object A must be updatedif ~get(A,'residflag'), A = residual(A); endif nargin ~= 2   error('Must specify paramlist')   end   paramindex = keyindex(get(A,'param'),Aparamkey);if paramindex == 0, error('Invalid parameter'), endR = columnnae(A,'resid');P = get(A,'param',Aparamkey);if isa(P,'scalarfield')   P = get(P,'val');   endP = P(1:prod(size(P)),1);Epsil = sqrt(eps)*max([abs(P),ones(size(P))],[],2).*sgn(P);PPlusEpsil = P + Epsil;% obsolete...% Epsil = 1e-4*mean(abs(P));% if Epsil < sqrt(eps), Epsil = sqrt(eps); end% PPlusEpsil = P*(1+Epsil);% PDiff = Epsil*P;% I = find( abs(PDiff) < Epsil );% PPlusEpsil(I) = P(I) + Epsil;% PDiff(I) = Epsil;n = 0;B = A;% Compute sensitivities of residuals to all keys in% cell array Akey; each (set of) column(s) corresponds% to each parameter in Akey.for j = 1:length(P)   B = set(B,'param',values( setval(get(A,'param'),PPlusEpsil(j),Aparamkey,j) ));   B = residual(B);   RPlusEpsil = columnnae(B,'resid');   n = n+1;   S(:,n) = ( RPlusEpsil - R )/Epsil(j);% obsolete S(:,n) = ( RPlusEpsil - R )/PDiff(paramindex);   end% -------------------- sgn.m -------------------- %function B = sgn(A)% return the sign of A if nonzero; return 1 if zeroI = find( A==0 );if ~isempty(I)   B(I,1) = 1;   endJ = find( A~=0 );if ~isempty(J)   B(J,1) = A(J,1)./abs(A(J,1));   end