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 

endfunction
