function [G, D, abstol] = computeObjSens(m, con, obj, opts)
verboseAll = max(opts.Verbose-1,0);

% Constants
nX = m.nX;
nVP = sum(sum(opts.UseParams));
nVX = sum(sum(opts.UseICs));
nV  = nVP + nVX;
nCon = length(con);
nObj = size(obj,2);

% Initialize variables
G = 0;
D = zeros(nV,1);
VXind = 1; % Stores the position in D where the first variable IC goes for each iCon
intOpts = opts;
if nargout >= 3; abstol = cell(nCon,1); end

if opts.Verbose; disp('Integrating sensitivities:'); end
for iCon = 1:nCon
    if verboseAll; tic; end
    
    % 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};
    intOpts.ObjWeights = opts.ObjWeights(iCon,:);
    intOpts.tGet = opts.tGet{iCon};
    
    % Integrate
    if opts.continuous(iCon) && opts.complex(iCon)
        sol = integrateObjSens(m, con(iCon), obj(iCon,:), intOpts);
    elseif opts.complex(iCon)
        sol = integrateSens(m, con(iCon), intOpts);
    elseif opts.continuous(iCon)
        sol = integrateObjSensDisc(m, con(iCon), obj(iCon,:), intOpts);
    else
        sol = integrateSensDisc(m, con(iCon), intOpts);
    end
    
    % *Compute G*
    % Extract continuous term
    if opts.continuous(iCon)
        contG = sol.y(nX+1,end);
    else
        contG = 0;
    end
    
    % Compute discrete term
    discG = 0;
    discreteTimes = [];
    for iObj = 1:nObj
        [iDiscG, temp] = obj(iCon,iObj).G([], sol, con(iCon).u, 1);
        discreteTimes = [discreteTimes; vec(temp)];
        discG = discG + opts.ObjWeights(iCon,iObj) * iDiscG;
    end

    % Remove repetitive discreteTimes
    discreteTimes = unique(discreteTimes);
    nDisc = numel(discreteTimes);

    % Add to cumulative goal value
    G = G + contG + discG;
    
    % *Compute D*
    % Extract continuous term
    if opts.continuous(iCon)
        dGdvStart = nX+1+nX*inV+1;
        dGdvEnd   = nX+1+nX*inV+inV;
        contD = sol.y(dGdvStart:dGdvEnd,end);
    else
        contD = sparse(nV,1);
    end
    
    % Compute discrete term
    dxdvStart = nX+opts.continuous(iCon)+1;
    dxdvEnd   = nX+opts.continuous(iCon)+nX*inV;
    dxdv = deval(sol, discreteTimes, dxdvStart:dxdvEnd); % xv_t
    discD = zeros(inV,nObj*nDisc); % v_ot
    for iObj = 1:nObj
        objDiscD = zeros(inV,1); % v_
        for iDisc = 1:nDisc
            objDiscD = objDiscD + vec(obj(iCon,iObj).dGdx(discreteTimes(iDisc), sol, con(iCon).u, 1) * reshape(dxdv(:,iDisc), nX, inV)); % v_ + (_x * x_v --> _v --> v_) --> v_
            temp = vec(obj(iCon,iObj).dGdp(discreteTimes(iDisc), sol, con(iCon).u, 1)); % p_ % partial dGdp(i)
            temp = [temp(opts.UseParams); sparse(inVX,1)]; % p_ --> v_
            objDiscD = objDiscD + temp; % v_
        end
        discD(:,(iObj-1)*nDisc + iDisc) = opts.ObjWeights(iCon,iObj) * objDiscD; % v_ as a row in v_ot
    end
    
    % Sorting by abs value to minimize numerical error
    [unused I] = sort(abs(discD),2);
    for i = 1:inV
        discD(i,:) = discD(i,I(i,:));
    end
    discD = sum(discD,2);
    
    % Sum discrete and continuous terms
    if opts.UseModelICs
        % All conditions have the same variables
        D = D + contD + discD;
    else
        % Rate parameters are the same
        D(1:nVP) = D(1:nVP) + contD(1:nVP) + discD(1:nVP);
        % IC parameters are different
        D(VXind:VXind+inVX-1) = contD(nVP+1:nVP+inVX) + discD(nVP+1:nVP+inVX);
        % Increment index of variable ICs
        VXind = VXind + inVX;
    end
    
    % Return the adaptive abstol, if requested
    if nargout >= 3
        abstol{iCon} = opts.RelTol * abstolObjSimple(m, con(iCon), obj(iCon,:), sol);
    end
    
    if verboseAll; fprintf('iCon = %d\t||dGdp|| = %g\tTime = %0.2f\n', iCon, norm(contD + discD), toc); end
end

if opts.Verbose;fprintf('Summary: ||dGdp|| = %g\n', norm(D));end
