function [varargout] = SimulateSensitivity(m, con, opts)
%SIMULATESENSITIVTY Integrate the sensitivities of every species with
%   respect to every parameter over all time. This function computes the
%   derivative on the species concentration with respect to the parameters
%   over all time.
%
%   Mathematically: dx/dp = Integral(df/dx * dx/dp + df/dp)
%   
%   [...] = SimulateSensitivity(m, con, opts)
%   
%   Inputs:
%       m    - The KroneckerBio model that will be simulated
%       con  - An array of experimental conditions under which the model
%              will be simulated
%       opts - Optional function 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 = []
%           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 
%           tGet        - Numerical vector for what time points will be 
%                         returned. Empty vector returns a solution that
%                         can be evalutated at all time points. Indicating
%                         tGet saves memory. Default = []
%           Verbose     - Numerical scalar for how much progress 
%                         information should be displayed. Default = 0
%
%   Outputs:
%       SimulateSensitivity(m, con, opts)
%        - Plots the sensitivities under each condition
%
%       sim = SimulateSensitivity(m, con, opts)
%        - An array of structures with each entry being the sensitivities
%          under one of the conditions.
%           sim.sol - The integrator solution to the sensitivies.
%           sim.t   - If opts.tGet is empty: A vector of timepoints chosen 
%                     by the ode solver.
%                     If opts.tGet is specified: This is opts.tGet.
%           sim.y   - If opts.tGet is empty: A function handle @(t, y) that
%                     evaluates the sensitivities at some particular
%                     timepoints (t) with respect to some particular output
%                     indexes (y).
%                     If opts.tGet is specified: A matrix m.nY by length of
%                     tGet. Each column is the sensitivities at the
%                     corresponding tGet.
%           sim.x   - If opts.tGet is empty: A function handle @(t, x) that
%                     evaluates the sensitivities at some particular
%                     timspoints (t) with respect to some particular
%                     species indexes (x).
%                     If opts.tGet is specified: A matrix m.nX by length of
%                     tGet. Each column is the sensitivties at the
%                     corresponding tGet.
%           sim.te  - Times at which each event occured
%           sim.ye  - Value of outputs at each event
%           sim.xe  - Value of species at each event
%           sim.ie  - Index of each event
%       
%       [t, y] = SimulateSensitivity(m, con, opts)
%        - sim.t, sim.y
%
%       [t, y, x] = SimulateSensitivity(m, con, opts)
%        - sim.t, sim.y, sim.x
%
%       [t, y, x, sol] = SimulateSensitivity(m, con, opts)
%        - sim.t, sim.y, sim.x, sim.sol
%       
%	Special:
%       This function can also be requested to return an empty array of
%       structures with the same fields as simulation. This may be
%       necessary to initialize an array that the user will later fill.
%
%       empty = SimulateSensitivity(m)
%           For integer m, create an empty simulation array m by 1
%
%       empty = SimulateSensitivity(m,n)
%           For integers m and n, create an empty simulation array m by n
%
%       empty = SimulateSensitivity([m,n,p...])
%           For integers array, create an empty simulation array size
%           [m,n,p...]

% (c) 2010 David R Hagen, Joshua F Apgar, Jared E Toettcher, & Bruce Tidor
% This work is released under the MIT license.

%% Work-up
% Clean up inputs
assert(nargout <= 4, 'KroneckerBio:SimulateSensitivity:FourOrFewerOutputs', 'SimulationGradient must have between 0 and 4 outputs.')
if nargin < 3
    opts = [];
    if nargin < 2
        con = [];
        if nargin < 1
            m = [];
        end
    end
end

% Special case: return empty structure array if inputs are numeric
if isnumeric(m)
    assert(nargout <= 1, 'KroneckerBio:SimulateSensitivity:OneOutputOnEmpty', 'SimulationGradient must have only one output when an empty simulation structure is requested.')
    temp = struct('sol',[],'t',[],'y',[],'x',[],'te',[],'ye',[],'xe',[],'ie',[]);
    if isempty(m)
        m = 1;
    end
    if ~isscalar(m)
        varargout{1}(prod(m),1) = temp;
        varargout{1} = reshape(varargout{1}, m);
        return
    end
    if isempty(con)
        con = 1;
    end
    varargout{1}(m, con) = temp;
    return
end

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

% Options
defaultOpts.UseParams   = 1:m.nP;
defaultOpts.UseICs      = [];
defaultOpts.UseModelICs = true;
defaultOpts.AbsTol      = NaN;
defaultOpts.RelTol      = NaN;
defaultOpts.tGet        = [];
defaultOpts.Verbose     = 0;

opts = mergestruct(defaultOpts, opts);

verbose = logical(opts.Verbose);
opts.Verbose = max(opts.Verbose-1,0);

% Constants
nX = m.nX;
nY = m.nY;
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);

nV = nVX + nVP;

% 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, false(nCon,1), nX, nCon, false, opts.UseParams, opts.UseICs, opts.UseModelICs);

%% Run integration for each experiment
sim(nCon) = struct('sol',[],'t',[],'y',[],'x',[],'te',[],'ye',[],'xe',[],'ie',[]);
intOpts = opts;

for iCon = 1:nCon
    % If opts.UseModelICs is false, the number of variables can change
    if opts.UseModelICs
        inVX = nVX;
        inV = nV;
    else
        inVX = sum(opts.UseICs(:,iCon));
        inV = nVP + inVX;
    end
	
    % Modify opts structure
    intOpts.AbsTol = opts.AbsTol{iCon};
    
    % Integrate dx/dp over time
    if verbose; fprintf(['Integrating sensitivities for ' con(iCon).name '...']); end
    if isempty(opts.tGet)
        sol = integrateSens(m, con(iCon), intOpts);
    else
        sol = integrateSensDisc(m, con(iCon), intOpts);
    end
    if verbose; fprintf('done.\n'); end
    
    % Store results
    if isempty(opts.tGet)
        % Return full solution
        sim(iCon).sol    = sol;
        sim(iCon).sol.c  = m.c;
        sim(iCon).t      = sol.x;
        sim(iCon).y      = @(varargin)evaluateOutputs(sol, varargin{:});
        sim(iCon).x      = @(varargin)evaluateSpecies(sol, varargin{:});
        sim(iCon).te     = sol.xe;
        if ~isempty(sol.ye)
            sim(iCon).ye = m.c*sol.ye;
        else
            sim(iCon).ye = [];
        end
        sim(iCon).xe     = sol.ye;
        sim(iCon).ie     = sol.ie;
    else
        % Return discrete solution
        sim(iCon).sol    = sol;
        sim(iCon).sol.c  = m.c;
        sim(iCon).t      = sol.x;
        sim(iCon).y      = reshape(m.c*reshape(sol.y, nX,(inV+1)*length(sol.x)), nY*(inV+1),length(sol.x));
        sim(iCon).x      = sol.y;
        sim(iCon).te     = sol.xe;
        if ~isempty(sol.ye)
            sim(iCon).ye = m.c*sol.ye;
        else
            sim(iCon).ye = [];
        end
        sim(iCon).xe     = sol.ye;
        sim(iCon).ie     = sol.ie;
    end
end

%% Work-down
switch (nargout)
    case 0
        % Save the hold state of the figure
        holdState = ishold;
        % Draw each result
        for iCon = 1:nCon
            plotExperiment(m, sim(iCon));
            hold on;
        end
        % Reset the hold state
        if ~holdState
            hold off;
        end
    case 1
        varargout{1} = sim;
    case 2 
        varargout{1} = sim.t;
        varargout{2} = sim.y;
    case 3
        varargout{1} = sim.t;
        varargout{2} = sim.y;
        varargout{3} = sim.x;
    case 4
        varargout{1} = sim.t;
        varargout{2} = sim.y;
        varargout{3} = sim.x;
        varargout{4} = sim.sol;
end

% End of function
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Evaluation functions %%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    function val = evaluateOutputs(sol, t, ind)
        nV = (size(sol.y, 1) - nX) / nX;
        nT = length(sol.x);
        
        if nargin < 3
            val = deval(sol, t, nX+1:nX+nX*nV); % xp_t
            val = reshape(val, nX,nV*nT); % x_pt
            val = m.c*val; % y_pt
            val = reshape(val, nY*nV,nT); % yp_t
        else
            if islogical(ind)
                % Convert to linear indexes
                ind = find(ind);
            end
            val = deval(sol, t, nX+1:nX+nX*nV); % xp_t
            val = reshape(val, nX,nV*nT); % x_pt
            val = m.c(ind,:)*val; % y_pt
            val = reshape(val, numel(ind)*nV,nT); % yp_t
        end
    end

    function val = evaluateSpecies(sol, t, ind)
        nV = (size(sol.y, 1) - nX) / nX;
        nT = length(sol.x);
        
        if nargin < 3
            val = deval(sol, t, nX+1:nX+nX*nV); % xp_t
        else
            if islogical(ind)
                % Convert to linear indexes
                ind = find(ind);
            end
            val = deval(sol, t, nX+1:nX+nX*nV); % xp_t
            val = reshape(val, nX,nV*nT); % x_pt
            val = val(ind,:); % x_pt chopped out rows
            val = reshape(val, numel(ind)*nV,nT); % xp_t
        end
    end

end