function H = sfd(x,grad,H,group,alpha,DiffMinChange,DiffMaxChange, ...
                  funfcn,varargin)
%SFD	Sparse Hessian via finite gradient differences
%
% H = sfd(x,grad,H,group,alpha,DiffMinChange,DiffMaxChange,funfcn,varargin)
% returns the sparse finite difference approximation H of the Hessian matrix 
% of the function 'funfcn' at the current point x. The vector group indicates 
% how to use sparse finite differencing: group(i) = j means that column i 
% belongs to group (or color) j. Each group (or color) corresponds to a finite 
% gradient difference.
%
% A non-empty input alpha overrides the default finite differencing 
% stepsize.
% 

%   Copyright 1990-2004 The MathWorks, Inc.
%   $Revision: 1.5.4.2 $  $Date: 2004/02/07 19:13:55 $

xcurr = x(:);  % Preserve x so we know what funfcn expects
scalealpha = false;
[m,n] = size(H); 
v = zeros(n,1);
ncol = max(group); 
if isempty(alpha)
    scalealpha = true;
    alpha = repmat(sqrt(eps),ncol,1);
end
H = spones(H); 
for k = 1:ncol
    d = (group == k);
    if scalealpha 
        xnrm = norm(xcurr(d));
        xnrm = max(xnrm,1);
        alpha(k) = alpha(k)*xnrm;
    end

    % Ensure magnitude of step-size lies within interval [DiffMinChange, DiffMaxChange]
    alpha(k) = sign(alpha(k))*min(max(abs(alpha(k)),DiffMinChange), ...
                                  DiffMaxChange);
    
    y = xcurr + alpha(k)*d;    
    
    % Make x conform to user-x
    x(:) = y;

    switch funfcn{1}
    case 'fun'
        error('optim:sfd:ShouldNotReachThis','should not reach this')
    case 'fungrad'
        [dummy,v(:)] = feval(funfcn{3},x,varargin{:});
    case 'fun_then_grad'
        v(:) = feval(funfcn{4},x,varargin{:});
    otherwise
        error('optim:sfd:UndefCallType','Undefined calltype in FMINUNC');
    end
    
    
    w = (v-grad)/alpha(k);
    cols = find(d); 
    A = sparse(m,n);
    A(:,cols) = H(:,cols);
    H(:,cols) = H(:,cols) - A(:,cols);
    [i,j,val] = find(A);
    [p,ind] = sort(i);
    val(ind) = w(p);
    A = sparse(i,j,full(val),m,n);
    H = H + A;
end
H = (H+H')/2;   % symmetricize