function [logL,dlogL] = likelihood_all_ErbB(tin,theta,C,I,f_extract)

f_plot = 0;

% Generate Time Vector

if(isempty(tin))
    tout = [1     10     20     30     40     50    100    200    300    400    500   1000   2000   3000   4000   5000  10000];
else
    tout = tin;
end

logLj = zeros(1,length(I));
dlogLj = zeros(length(theta),length(I));

options.atol = 1e-8;
options.rtol = 1e-8;
options.maxsteps = 1e4;
options.linsol = 9;
options.sensi_meth = 'adjoint';
if(nargout>1)
    options.sensi = 1;
else
    options.sensi = 0;
end

np = size(theta,1);

nderiv = nargout;

D = amidata(1,1,0,0,size(C,1));
D.Y = 1;

if(f_extract>1)
    nderiv = 2;
    options.sensi = 1;
end

Iret = I;

try
    p = gcp('nocreate');
catch
    p = [];
end
if(isempty(p))
    for icl = 1:length(I)
        if(f_extract)
            I_icl = likelihood_cell_line(tout,theta,C,I,D,options,f_plot,f_extract,nderiv,icl);
            Iret{icl} = I_icl{icl};
        else
            [logL_icl,dlogL_icl] = likelihood_cell_line(tout,theta,C,I,D,options,f_plot,f_extract,nderiv,icl);
            logLj(:,icl) = logL_icl;
            if(nderiv>1)
                dlogLj(:,icl) = dlogL_icl;
            end
        end
    end
else
    parfor icl = 1:length(I)
        if(f_extract)
            I_icl = likelihood_cell_line(tout,theta,C,I,D,options,f_plot,f_extract,nderiv,icl);
            Iret{icl} = I_icl{icl};
        else
            [logL_icl,dlogL_icl] = likelihood_cell_line(tout,theta,C,I,D,options,f_plot,f_extract,nderiv,icl);
            logLj(:,icl) = logL_icl;
            if(nderiv>1)
                dlogLj(:,icl) = dlogL_icl;
            end
        end
    end
end

if(f_extract)
    logL = Iret;
else
    logL = sum(logLj,2);
    if(nderiv>1)
        dlogL = sum(dlogLj,2);
    end
end
end



function varargout = likelihood_cell_line(tout,theta,C,I,D,options,f_plot,f_extract,nderiv,icl)

logLj = zeros(1,1);
dlogLj = zeros(length(theta),1);

C_ixp = C(:,I{icl}.ref);
sol_ref = simulate_model_ErbB(tout,theta,C_ixp,D,options);
if(sol_ref.status<0)
    error('failed to integrate reference');
end
nx = size(sol_ref.x,2);
nc = size(C,1);
while(max(abs(sol_ref.xdot./(sol_ref.x+1e-6)))>1e-6)
    C_ixp((nc-nx+1):end) = sol_ref.x;
    sol_ref = simulate_model_ErbB(tout,theta,C_ixp,D,options);
    if(sol_ref.status<0)
        error('failed to integrate reference');
    end
end

for ixp = 1:length(I{icl}.exp)
    Y = I{icl}.D(ixp).Y;
    Sigma_Y = I{icl}.D(ixp).Sigma_Y;
    C_ixp = C(:,I{icl}.exp(ixp));
    C_ixp((nc-nx+1):end) = sol_ref.x;
    sol_drug = simulate_model_ErbB(tout,theta,C_ixp,I{icl}.D(ixp),options);
    if(sol_drug.status<0)
        error('failed to integrate drug');
    end
    while(max(abs(sol_drug.xdot./(sol_drug.x+1e-6)))>1e-6)
        C_ixp((nc-nx+1):end) = sol_drug.x;
        sol_drug = simulate_model_ErbB(tout,theta,C_ixp,I{icl}.D(ixp),options);
        if(sol_drug.status<0)
            error('failed to integrate drug');
        end
    end
    inhib = (sol_drug.llh/sol_ref.llh);
    if(f_extract)
        I{icl}.ref_ss = sol_ref.x;
        I{icl}.drug_ss(:,ixp) = sol_drug.x;
        I{icl}.sim(ixp) = inhib;
    end
    logLj = logLj - 0.5*( ( inhib - Y ) / Sigma_Y )^2;
    if(nderiv>1)
        dinhib = (sol_drug.sllh.*sol_ref.llh - sol_drug.llh.*sol_ref.sllh)./(sol_ref.llh.^2);
        dlogLj = dlogLj - ( ( inhib - Y )*( dinhib ) / Sigma_Y^2 );
        if(f_extract)
            if(f_extract>1)
                I{icl}.ssim(:,ixp) = dinhib;
            end

        end
    end
    
end

if(f_extract)
    varargout{1} = I;
else
    varargout{1} = logLj;
    varargout{2} = dlogLj;
end

end
