function df = d(func,k,varargin)
% d  [Not a public function] Compute numerical derivatives of non-analytical or user-defined functions.
%
% Backend IRIS function.
% No help provided.

% -IRIS Toolbox.
% -Copyright (c) 2007-2012 Jaromir Benes.

if ~isempty(varargin) ...
        && (iscellstr(varargin{end}) || isequal(varargin{end},Inf))
    userdifflist = varargin{end};
    varargin(end) = [];
else
    userdifflist = {};
end

%**************************************************************************

% Try to get a user-supplied first derivative by calling the
% function itself with a 'diff' argument. For example,
%
%     func(x,y,z,'diff',3)
%
% expects the first derivative of the function w.r.t. the 3-rd input
% argument.

nd = length(varargin{k});
if any(strcmp(char(func),userdifflist)) || isequal(userdifflist,Inf)
    status = warning();
    warning('off'); %#ok<WNOFF>
    try
        % User-supplied derivatives.
        df = feval(func,varargin{:},'diff',k);
        if ~isnumeric(df) || length(df) ~= nd
            df = NaN;
        end
    catch %#ok<CTCH>
        df = NaN;
    end
    warning(status);
    if isfinite(df)
        return
    end
end

% Compute the derivative numerically.
if length(k) == 1
    % First derivative.
    df = xxdiff(func,k,varargin{:});
    return
elseif length(k) == 2
    % Second derivative; these are needed in optimal policy models with
    % user-supplied functions.
    y0 = varargin{k(2)};
    hy = abs(eps()^(1/3.5))*max([y0,1]);
    yp = y0 + hy;
    ym = y0 - hy;
    varargin{k(2)} = yp;
    fp = xxdiff(func,k(1),varargin{:});
    varargin{k(2)} = ym;
    fm = xxdiff(func,k(1),varargin{:});
    df = (fp - fm) / (yp - ym);
end

end

% Subfunctions.

%**************************************************************************
function df = xxdiff(func,k,varargin)

epsilon = eps()^(1/3.5);
x0 = varargin{k};
hx = abs(epsilon*max(x0,1));
xp = x0 + hx;
xm = x0 - hx;
varargin{k} = xp;
fp = feval(func,varargin{:});
varargin{k} = xm;
fm = feval(func,varargin{:});
df = (fp - fm) ./ (xp - xm);

end
% xxdiff().