function p = ObjectiveProbability(m, con, obj, opts)
% p = ObjectiveProbability(m, con, obj, opts)

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

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

% Options
defaultOpts.RelTol      = NaN; % 1e-6
defaultOpts.AbsTol      = NaN; % in fixAbsTol
defaultOpts.UseModelICs = true;
defaultOpts.UseAdjoint  = false;

defaultOpts.Verbose     = 0;

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

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

% Standardize structures
con = pastestruct(Uzero(m), con);
obj = pastestruct(Gzero(m), obj);

%% Fetch 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, 1, opts.continuous, nx, nCon);

%% Probability
% Initialize variables
p = 1;

for iCon = 1:nCon
    if verbose; disp(['Integrating system for ' con(iCon).name '...']); end
    % Modify opts structure
    optsSub            = opts;
    optsSub.AbsTol     = opts.AbsTol{iCon};
    optsSub.tGet       = opts.tGet{iCon};
    
    % Integrate
    if opts.continuous(iCon) && opts.complex(iCon)
        sol = integrateObj(m, con(iCon), obj(iCon,:), optsSub);
    elseif opts.complex(iCon)
        sol = integrateSys(m, con(iCon), optsSub);
    elseif opts.continuous(iCon)
        sol = integrateObjDisc(m, con(iCon), obj(iCon,:), optsSub);
    else
        sol = integrateSysDisc(m, con(iCon), optsSub);
    end
    
    % Probability
    for iObj = 1:nObj
        p = p * obj(iCon,iObj).p(sol, con(iCon).u);
    end
end
