function [D All] = ObjectiveGradient(m, con, obj, opts)
%OBJECTIVEGRADIENT Evaluate the gradient of an objective function
%
%   Mathematically: dG/dp = sum(dG/dx(t) * dx/dp, t)
%
%   [D All] = ObjectiveGradient(m, con, obj, opts)
%
%   Inputs:
%       m    - The KroneckerBio model that will be simulated
%       con  - An vector of experimental conditions under which the model
%              will be simulated
%       obj  - A matrix length of con by any number of objective functions
%              that will be evaluated
%       opts - Scalar structure of options
%           UseModelICs - Logical scalar indicating that the model's
%                         initial conditions should be used instead of
%                         those of the experimental conditions.
%                         Default = true
%           UseParams   - Vector of indexes indicating the rate constants
%                         whose sensitivities will be considered.
%                         Default = 1:m.nP
%           UseICs      - Vector of indexes indicating the initial
%                         concentrations whose sensitivities will be
%                         considered. Default = []
%           Normalized  - Scalar logical indicating that the gradient
%                         should be normalized with respect to the
%                         parameter values. Default = true
%           UseAdjoint  - Scalar logical indicating that the sensitivites 
%                         should be calculated via the adjoint method
%                         rather than the forward method. Default = true
%           RelTol      - Numerical scalar for the relative tolerance of 
%                         the integration. Default = 1e-6
%           AbsTol      - Cell vector of vectors, numerical scalar, or 
%                         numerical vector for the absolute tolerance of
%                         the integration. Default = 1e-9 
%           ObjWeights  - A matrix of size(obj) indicating a post-
%                         evaluation weight on each objective function in
%                         terms of how much it will contribute to the final
%                         value. Default = ones(size(obj))
%           Verbose     - Numerical scalar for how much progress 
%                         information should be displayed. Default = 0
%
%   Outputs:
%       D   - The sum of all objective function gradients
%       All - A matrix size(obj) for the values for each objective function

% (c) 2010 David R Hagen & Bruce Tidor
% This work is released under the MIT license.

%% Work-up
% Clean up inputs
if nargin < 4
    opts = [];
end

assert(isscalar(m), 'KroneckerBio:ObjectiveGradient:MoreThanOneModel', 'The model structure must be scalar')

% Options
defaultOpts.UseModelICs             = true;
defaultOpts.UseParams               = 1:m.nP;
defaultOpts.UseICs                  = [];
defaultOpts.Normalized              = true;
defaultOpts.RelTol                  = NaN;
defaultOpts.AbsTol                  = NaN;
defaultOpts.ObjWeights              = ones(size(obj));
defaultOpts.UseAdjoint              = true;
defaultOpts.Verbose                 = 0;

opts = mergestruct(defaultOpts, opts);

% Constants
nX = m.nX;
nP = m.nP;
nCon = numel(con);

% Ensure UseRates is column vector of logical indexes
[opts.UseParams, nVP] = fixUseParams(opts.UseParams, nP);

% Ensure UseICs is a matrix of linear indexes
[opts.UseICs, nVX] = fixUseICs(opts.UseICs, opts.UseModelICs, nX, nCon);

% Add missing fields to structure
con = pastestruct(Uzero(m), con);
obj = pastestruct(Gzero(m), obj);

% Fix integration type
[opts.continuous, opts.complex, opts.tGet] = fixIntegrationType(obj);

%% Tolerances
% RelTol
if isempty(opts.RelTol) || isnan(opts.RelTol)
    opts.RelTol = 1e-6;
end

% Fix AbsTol to be a cell array of vectors appropriate to the problem
opts.AbsTol = fixAbsTol(opts.AbsTol, 2, opts.continuous, nX, nCon, opts.UseAdjoint, opts.UseParams, opts.UseICs, opts.UseModelICs);

%% Run appropriate integration
if nargout <= 1
    if opts.UseAdjoint
        [unused D] = integrateObjSensAdj(m, con, obj, opts);
    else
        [unused D] = computeObjSens(m, con, obj, opts);
    end
else
    if opts.UseAdjoint
        [D All] = integrateObjSensAdjAll(m, con, obj, opts);
    else
        [D All] = computeObjSensAll(m, con, obj, opts);
    end
end

%% Normalization
if opts.Normalized
    % Extract parameter set
    pAll = m.p;
    
    if opts.UseModelICs
        icAll = m.ic;
    else
        icAll = zeros(nX, nCon);
        for i = 1:nCon
            icAll(:,i) = con(i).ic;
        end
    end
    
    v = [pAll(opts.UseParams); vec(icAll(opts.UseICs))];
    
    % Normalize
    D = D .* v;
    if nargout >= 2
        for iCon = 1:nCon
            for iObj = 1:nObj
                All{iCon,iObj} = All{iCon,iObj} .* v;
            end
        end
    end
end