function [F All] = ObjectiveInformation(m, con, obj, opts, dxdpSol)
%OBJECTIVEINFORMATION Compute the fisher information matrix of an objective
%   function.
%
%   Given a model and experimental conditions, this function computes the
%   fisher information matrix provided by the objective function, given the
%   model and conditions and evaluated at p.
%
%   Mathematically: F = dy/dp' * V^-1 * dy/dp
%
%   [...] = ObjectiveInformation(m, con, obj, opts, dxdpSol)
%
%   Inputs:
%       m    - The KroneckerBio model that will be simulated
%       con  - A structure vector of the experimental conditions under
%              which the hessian will be evaluated
%       obj  - A structure array of the objective functions under which the
%              hessian will be evaluated. Each row of obj is matched to the
%              corresponding entry in con.
%       opts - Optional function options
%           .UseParams   - Vector of indexes indicating the rate constants
%                          whose sensitivities will be considered
%           .UseICs      - Vector of indexes indicating the initial
%                          concentrations whose sensitivities will be
%                          considered
%           .UseModelICs - Boolean that determines whether to use the
%                          initial concentrations of the model or the
%                          conditions. Default = true
%           .Normalized  - Logical that determines if the simple
%                          information or the normalized information will
%                          be computed. The normalized information is
%                          normalized with respect to the values of the
%                          parameters. Default = true
%           .Verbose     - Print progress to command window
%       dxdpSol   - A structure vector containing the solution to the model
%                   sensitivities under each condition. Optional, but
%                   speeds up the calculation.
%
%   Outputs:
%       F = ObjectiveInformation(m, con, obj, ...)
%           F - A matrix
%
%       [F, All] = ObjectiveInformation(m, con, obj, ...)
%           All - The fisher information matrix (FIM) is the sum of all
%                 fisher information matrices assuming there is no
%                 covariance between errors in seperate experiments. All
%                 provides the individual FIMs for each experiment.
%
%   Additional info:
%   - The experimental condition vector can also be a cell vector
%   - The objective function array can also be a cell array. Empty entries 
%   in the cell array and entries in the structure array with empty
%   values are ignored. This way, conditions can have different numbers of
%   objective functions associated with them.
%   - The optional solutions, dxdpSol, can also be cell vectors. Empty
%   entries in the cell vector and entries in the structure vector with
%   empty values are ignored.

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

% TODO: make the discreteTimes more general.

%% Work-up
% Clean up inputs
assert(nargin >= 3, 'KroneckerBio:ObjectiveInformation:AtLeastThreeInputs', 'ObjectiveExpectedHessian requires at least 3 input arguments.')
if nargin < 5
    dxdpSol = [];
    if nargin < 4
        opts = [];
    end
end

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

% Options
defaultOpts.UseParams   = 1:m.nP;
defaultOpts.UseICs      = [];
defaultOpts.UseModelICs = true;
defaultOpts.Normalized  = true;
defaultOpts.Verbose     = 0;

opts = mergestruct(defaultOpts, opts);

% Constants
nX = m.nX;
nP = m.nP;
nCon = length(con);
nObj = size(obj, 2);

% 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);

simOpts = opts;

nV = nVX + nVP;

% Make conditions consistent as cell vectors
if isstruct(con)
    con = num2cell(con);
end

% Make objectives consistent as cell arrays
if isstruct(obj)
    obj = num2cell(obj);
end

% Make solutions consistent as cell vectors
if isstruct(dxdpSol)
    dxdpSol = num2cell(dxdpSol);
end

%% Loop through conditions
F = zeros(nV,nV);
VXind = nVP+1; % Stores the position in F where the first variable IC goes for each iCon

% Initialize All array if requested
if nargout >= 2
    All = cell(nCon,nObj);
end

for iCon = 1:nCon
    % Compute IC size this iteration
    if ~opts.UseModelICs
        inVX = sum(opts.UseICs(:,iCon));
    end
    
    % Collect the discrete timepoints to capture
    simOpts.tGet = [];
    for iObj = 1:nObj
        if ~isempty(obj{iCon,iObj}) 
            if ~isfield(obj{iCon,iObj}, 'discreteTimes') || isempty(obj{iCon,iObj}.discreteTimes)
                simOpts.tGet = [];
                break
            else
                simOpts.tGet = unique([simOpts.tGet, vec(obj{iCon,iObj}.discreteTimes)']);
            end
        end
    end
    
    % Sensitivity integration if not provided
    if isempty(dxdpSol) || isempty(dxdpSol{iCon})
        if ~opts.UseModelICs
            simOpts.UseICs = opts.UseICs(:,iCon);
        end
        
        sim = SimulateSensitivity(m, con{iCon}, simOpts);
        dxdpSoliCon = sim.sol;
        dxdpSoliCon.C = m.c;
    else
        dxdpSoliCon = dxdpSol{iCon};
    end
    
    % Sum all FIMs as computed by each objective function
    if opts.Normalized
        if opts.Verbose; fprintf(['Computing fisher information matrix (normalized) for ' con{iCon}.name '...']); end
        % Loop through each objective function in the current row
        for iObj = 1:nObj
            % Ignore empty structures
            if ~isempty(obj{iCon,iObj}) && ~isempty(obj{iCon,iObj}.Fn)
                if opts.UseModelICs
                    p = [m.p(opts.UseParams); m.ic(opts.UseICs)]; % model initial conditions
                    curF = obj{iCon,iObj}.Fn(dxdpSoliCon, p);
                    F = F + curF;
                else
                    p = [m.p(opts.UseParams); con{iCon}.ic(opts.UseICs(:,iCon))]; % con initial conditions
                    curF = obj{iCon,iObj}.Fn(dxdpSoliCon, p);
                    % Add to correct place in F
                    linInd = [1:nVP, VXind:VXind+inVX-1]; % linear indexes of parameters
                    F(linInd,linInd) = F(linInd,linInd) + curF;
                end
                
                % Store FIM if requested
                if nargout >= 2
                    All{iCon, iObj} = curF;
                end
            end
        end
    else
        if opts.Verbose; fprintf(['Computing fisher information matrix (non-normalized) for ' con{iCon}.name '...']); end
        % Loop through each objective function in the current row
        for iObj = 1:nObj
            % Ignore empty structures
            if ~isempty(obj{iCon,iObj}) && ~isempty(obj{iCon,iObj}.F)
                if opts.UseModelICs
                    curF = obj{iCon,iObj}.F(dxdpSoliCon);
                    F = F + curF;
                else
                    curF = obj{iCon,iObj}.F(dxdpSoliCon);
                    % Add to correct place in F
                    linInd = [1:nVP, VXind:VXind+inVX-1]; % linear indexes of parameters
                    F(linInd,linInd) = F(linInd,linInd) + curF;
                end
                
                % Store FIM if requested
                if nargout >= 2
                    All{iCon, iObj} = curF;
                end
            end
        end
    end
    
    % Update condition IC position
    if ~opts.UseModelICs
        VXind = VXind + inVX;
    end
    
    if opts.Verbose; fprintf('done.\n'); end
end