function [gradf,cJac,NEWLAMBDA,OLDLAMBDA,s] = finitedifferences(xCurrent,...
             xOriginalShape,funfcn,confcn,lb,ub,fCurrent,cCurrent,XDATA,YDATA,...
             DiffMinChange,DiffMaxChange,typicalx,finDiffType,variables,...
             LAMBDA,NEWLAMBDA,OLDLAMBDA,POINT,FLAG,s,isFseminf,startnlineq,varargin)
%FINITEDIFFERENCES computes finite-difference derivatives.
%
% This helper function computes finite-difference derivatives of the objective 
% and constraint functions.
%
%  [gradf,cJac,NEWLAMBDA,OLDLAMBDA,s] = FINITEDIFFERENCES(xCurrent, ... 
%                  xOriginalShape,funfcn,confcn,lb,ub,fCurrent,cCurrent, ...
%                  XDATA,YDATA,DiffMinChange,DiffMaxChange,typicalx,finDiffType, ...
%                  variables,LAMBDA,NEWLAMBDA,OLDLAMBDA,POINT,FLAG,s, ...
%                  varargin)
% computes the finite-difference gradients of the objective and
% constraint functions.
%
%  gradf = FINITEDIFFERENCES(xCurrent,xOriginalShape,funfcn,[],lb,ub,fCurrent,...
%                  [],YDATA,DiffMinChange,DiffMaxChange,typicalx,finDiffType,...
%                  variables,[],[],[],[],[],[],varargin)
% computes the finite-difference gradients of the objective function.
%
%
% INPUT:
% xCurrent              Point where gradient is desired
% xOriginalShape        Shape of the vector of variables supplied by the user
%                       (The value of xOriginalShape is NOT used)
% funfcn, confcn        Cell arrays containing info about objective and
%                       constraints, respectively. The objective (constraint) 
%                       derivatives are computed if and only if funfcn 
%                       (confcn) is nonempty.
%                       
% lb, ub                Lower and upper bounds
% fCurrent, cCurrent    Values at xCurrent of the function and the constraints 
%                       to be differentiated. Note that fCurrent can be a scalar 
%                       or a (row or column) vector. 
%
% XDATA, YDATA          Data passed from lsqcurvefit.
%
% DiffMinChange, 
% DiffMaxChange         Minimum and maximum values of perturbation of xCurrent 
% finDiffType           Type of finite difference desired (only forward 
%                       differences implemented so far)
% variables             Variables w.r.t which we want to differentiate. Possible 
%                       values are 'all' or an integer between 1 and the
%                       number of variables.
%
% LAMBDA,NEWLAMBDA,
% OLDLAMBDA,POINT,
% FLAG,s                Parameters for semi-infinite constraints
%
% isFseminf             True if caller is fseminf, false otherwise
%
% startnlineq           Another parameter used fseminf's semicon
%
% varargin              Problem-dependent parameters passed to the objective and 
%                       constraint functions
%
% OUTPUT:
% gradf                 If fCurrent is a scalar, gradf is the finite-difference 
%                       gradient of the objective; if fCurrent is a vector,
%                       gradf is the finite-difference Jacobian  
% cJac                  Finite-difference Jacobian of the constraints
% NEWLAMBDA,
% OLDLAMBDA,s           Parameters for semi-infinite constraints

%   Copyright 1990-2007 The MathWorks, Inc.
%   $Revision: 1.1.4.9 $  $Date: 2007/12/10 22:34:07 $

% For vector-valued functions in funfcn, we make the function values 
% (both at xCurrent and those at the perturbed points) column vectors 
% to ensure that the given fCurrent and the computed fplus will have the 
% same shape so that fplus - fCurrent be well defined.  
fCurrent = fCurrent(:);
numberOfFunctions = numel(fCurrent);
numberOfVariables = numel(xCurrent); 
functionIsScalar = (numberOfFunctions == 1);

% nonEmptyLowerBounds = true if lb is not empty, false if it's empty;
% analogoulsy for nonEmptyUpperBound
nonEmptyLowerBounds = ~isempty(lb);
nonEmptyUpperBounds = ~isempty(ub);

% Make sure xCurrent and typicalx are column vectors so that the 
% operation max(abs(xCurrent),abs(typicalx)) won't error
xCurrent = xCurrent(:); typicalx = typicalx(:);
% Value of stepsize suggested in Trust Region Methods, Conn-Gould-Toint, section 8.4.3
CHG = sqrt(eps)*sign(xCurrent).*max(abs(xCurrent),abs(typicalx));
%
% Make sure step size lies within DiffminChange and DiffMaxChange
%
CHG = sign(CHG+eps).*min(max(abs(CHG),DiffMinChange),DiffMaxChange);
len_cCurrent = length(cCurrent); % For semi-infinite

if nargout < 3
   NEWLAMBDA=[]; OLDLAMBDA=[]; s=[];
end
if nargout > 1
      cJac = zeros(len_cCurrent,numberOfVariables);  
else
      cJac = [];
end
% allVariables = true/false if finite-differencing wrt to all/one variables
allVariables = false;
if ischar(variables)
   if strcmp(variables,'all')
      variables = 1:numberOfVariables;
      allVariables = true;
   else
      error('optimlib:finitedifferences:InvalidVariables', ...
            'Unknown value of input ''variables''.')
   end
end

% Preallocate gradf for speed 
if ~isempty(funfcn)
    if functionIsScalar
        gradf = zeros(numberOfVariables,1);
    elseif allVariables % vector-function and gradf estimates full Jacobian
        gradf = zeros(numberOfFunctions,numberOfVariables);
    else % vector-function and gradf estimates one column of Jacobian
        gradf = zeros(numberOfFunctions,1);
    end
else
    gradf = [];
end

% Do this switch outside of loop for speed
if isFseminf
   vararginAdjusted = varargin(3:end);
else
   vararginAdjusted = varargin;     
end

% If lsqcurvefit, then add XDATA to objective's input list.
% xargin{1} will be updated right before each evaluation 
if ~isempty(XDATA)
    xargin = {xOriginalShape,XDATA};
else
    xargin = {xOriginalShape};
end

for gcnt=variables
   if gcnt == numberOfVariables, 
      FLAG = -1; 
   end
   temp = xCurrent(gcnt);
   xCurrent(gcnt)= temp + CHG(gcnt);
         
   if (nonEmptyLowerBounds && isfinite(lb(gcnt))) || (nonEmptyUpperBounds && isfinite(ub(gcnt)))
      % Enforce bounds while finite-differencing.
      % Need lb(gcnt) ~= ub(gcnt), and lb(gcnt) <= temp <= ub(gcnt) to enforce bounds.
      % (If the last qpsub problem was 'infeasible', the bounds could be currently violated.)
      if (lb(gcnt) ~= ub(gcnt)) && (temp >= lb(gcnt)) && (temp <= ub(gcnt)) 
          if  ((xCurrent(gcnt) > ub(gcnt)) || (xCurrent(gcnt) < lb(gcnt))) % outside bound ?
              CHG(gcnt) = -CHG(gcnt);
              xCurrent(gcnt)= temp + CHG(gcnt);
              if (xCurrent(gcnt) > ub(gcnt)) || (xCurrent(gcnt) < lb(gcnt)) % outside other bound ?
                  [newchg,indsign] = max([temp-lb(gcnt), ub(gcnt)-temp]);  % largest distance to bound
                  if newchg >= DiffMinChange
                      CHG(gcnt) = ((-1)^indsign)*newchg;  % make sure sign is correct
                      xCurrent(gcnt)= temp + CHG(gcnt);
                      warning('optimlib:finitedifferences:StepReduced', ...
                             ['Derivative finite-differencing step was artificially reduced to be within\n', ...
                              'bound constraints. This may adversely affect convergence. Increasing distance between\n', ...
                              'bound constraints, in dimension %d to be at least %0.5g may improve results.'], ...
                              gcnt,abs(2*CHG(gcnt)))
                  else
                      error('optimlib:finitedifferences:DistanceTooSmall', ...
                          ['Distance between lower and upper bounds, in dimension %d is too small to compute\n', ...
                           'finite-difference approximation of derivative. Increase distance between these\n', ...
                           'bounds to be at least %0.5g.'],gcnt,2*DiffMinChange)
                  end          
              end
          end
      end
   end % of 'if isfinite(lb(gcnt)) || isfinite(ub(gcnt))'
   
   xOriginalShape(:) = xCurrent;
   xargin{1} = xOriginalShape; % update x in list of input arguments to objective
   if ~isempty(funfcn) % Objective gradient required
       % The length of varargin (depending on the caller being fseminf or not)
       % was calculated outside of the loop for speed.
       fplus = feval(funfcn{3},xargin{:},vararginAdjusted{:});
   else
       fplus = [];
   end
   % YDATA: Only used by lsqcurvefit, which has no nonlinear constraints
   % (the only type of constraints we do finite differences on: bounds 
   % and linear constraints do not require finite differences) and thus 
   % no needed after evaluation of constraints
   if ~isempty(YDATA)    
      fplus = fplus - YDATA;
   end
   % Make sure it's in column form
   fplus = fplus(:);
   if ~isempty(funfcn)
       if functionIsScalar
           gradf(gcnt,1) =  (fplus-fCurrent)/CHG(gcnt);
       elseif allVariables % vector-function and gradf estimates full Jacobian
           gradf(:,gcnt) = (fplus-fCurrent)/CHG(gcnt);
       else % vector-function and gradf estimates only one column of Jacobian
           gradf = (fplus-fCurrent)/CHG(gcnt);
       end
   end
   
   if ~isempty(cJac) % Constraint gradient required
         if isFseminf 
            [ctmp,ceqtmp,NPOINT,NEWLAMBDA,OLDLAMBDA,s] = ...
               semicon(xOriginalShape,LAMBDA,NEWLAMBDA,OLDLAMBDA,POINT,FLAG,s,startnlineq,varargin{:});            
         else
            [ctmp,ceqtmp] = feval(confcn{3},xOriginalShape,varargin{:});
         end
         cplus = [ceqtmp(:); ctmp(:)];

      % Next line used for problems with varying number of constraints
      if isFseminf && len_cCurrent~=length(cplus)
         cplus=v2sort(cCurrent,cplus); 
      end      
      if ~isempty(cplus)
         cJac(:,gcnt) = (cplus - cCurrent)/CHG(gcnt); 
      end           
   end
    xCurrent(gcnt) = temp;
end % for 








