function [max_errors, varargout] = var_checkgrad(fn, epsilon, varargin)
%VAR_CHECKGRAD numerically check function returning derivatives wrt its inputs
%
% [max_errors, varargout] = var_checkgrad(fn, epsilon, varargin)
%
% fn must be a function that will take all the varargin arguments (and no
% others). It should first return a scalar objective function value, followed
% by the partial derivatives for all of its input arguments.
%
% Each component of each member of varargin is perturbed by epsilon to
% numerically approximate the partial derivatives. The fractional max(abs(.))
% errors for each element of varargin is reported.
%
% The numerical derivatives are returned as a variable number of args if
% requested.

% Iain Murray, November 2007, July 2008, Feb 2015

assert(isscalar(epsilon), 'epsilon tolerance should be a scalar');

args = varargin;
grads = cell(size(args));
[cost, grads{:}] = fn(args{:});

for i = 1:length(args)
    assert(isequal(size(grads{i}), size(args{i})), ...
        sprintf('%dth gradient is wrong size', i));
end

max_errors = zeros(size(args));
diffs = grads;

for i = 1:length(args)
    for c = 1:prod(size(args{i}))
        args{i}(c) = args{i}(c) + epsilon/2;
        cost_f = fn(args{:});
        args{i}(c) = args{i}(c) - epsilon;
        cost_b = fn(args{:});
        diffs{i}(c) = (cost_f - cost_b)/epsilon;
        args{i}(c) = varargin{i}(c); % undo component tweak
    end
    max_errors(i) = max(abs(grads{i}(:)-diffs{i}(:))./grads{i}(:));
end

if nargout > 1
    varargout = diffs(1:(nargout-1));
end