saGradientMethod(ds,varargin)

function [gradf, ds] = saGradientMethod(ds,varargin)
% [gradf, ds] = saGradientMethod(ds,paramNames,responseNames)

%   Copyright 2007-2009 The MathWorks, Inc.
paramNames = {};
responseNames = {};
DiffMinChange = 1e-8;
scaled = 1;
ScaledMinChange = false;

% assign inputs
if nargin > 1, paramNames = varargin{1}; end
if nargin > 2, responseNames = varargin{2}; end
if nargin > 3, scaled = varargin{3}; ScaledMinChange = true; end

%% Pull out variables for Simulink to update
idx = ismember(ds.Properties.UserData, 'slParam');
Vars  = ds.Properties.VarNames(idx);

%% Pull out current value and bounds
if isempty(paramNames)
    paramNames = Vars;
end

xC = double(ds('Initial',paramNames));
lb = double(ds('lowerBound',paramNames));
ub = double(ds('upperBound',paramNames));

%% Run Simulink Model for current value
ds('Initial',:) = runModel(ds('Initial',:));
if isempty(responseNames)
    idx = ismember(ds.Properties.UserData, 'Response');
    responseNames  = ds.Properties.VarNames(idx);
end
noResponses = length(responseNames);
    
%% Calculate dX fore each variable in ds
TypicalX = ones(size(xC));
deltaX = (eps^(1/3))*max(abs(xC),abs(TypicalX));

for dim = 1:length(xC)
    if ScaledMinChange
        DiffMinChange = scaled*xC(dim);
    end
    if deltaX(dim) < DiffMinChange
        deltaX(dim) = DiffMinChange;
    end
    [deltaX(dim),formulaType] = cntrlFinDiffInsideBnds(xC(dim),lb(dim),...
        ub(dim),deltaX(dim),dim,DiffMinChange);
    
    xCElementOrig = xC(dim); % save component of current point
    
    % Perturb component of current point
    if formulaType == 0      % central difference formula
        xPert1Element = xCElementOrig - deltaX(dim);
        xPert2Element = xCElementOrig + deltaX(dim);
    elseif formulaType == -1 % double left step difference formula
        xPert1Element = xCElementOrig - 2*deltaX(dim);
        xPert2Element = xCElementOrig - deltaX(dim);
    else                     % double right step difference formula
        xPert1Element = xCElementOrig + deltaX(dim);
        xPert2Element = xCElementOrig + 2*deltaX(dim);
    end
    
    % Update dataset array with points to run
    xPert = ds(1,paramNames);
    xPert{1,paramNames{dim}} = xPert1Element;
    ds(['deltaX',num2str(dim),'1'],paramNames) = xPert;
    xPert{1,paramNames{dim}} = xPert2Element;
    ds(['deltaX',num2str(dim),'2'],paramNames) = xPert;
    
    % run Simulink Model
    ds(end-1:end,:) = runModel(ds(end-1:end,:));
           
    % calculate the gradient
    for i = 1:noResponses
        gradf(dim,1) = twoStepFinDiffFormulas(formulaType,deltaX(dim), ...
                        ds.(responseNames{i})('Initial'),...
                        ds.(responseNames{i})(end-1),...
                        ds.(responseNames{i})(end));
    end
end