%X(1) --> T   = Healthy T-cells
%X(2) --> T_I = Infected T-cells
%X(3) --> T_L = Latently-infected T-cells
%X(4) --> M   = Healthy macrophages
%X(5) --> M_I = Infected macrophages
%X(6) --> C   = Cytotoxic T-lymphocytes population
%X(7) --> V   = HIV population

%Initialize algorithm parameters
N = 1000;         %Number of samples
q = zeros(N,1);   %Output of interest (T-cells after 3000 days)
A = zeros(N,27);  %Matrix of random parameter samples

%Set upper and lower bounds for paramters
xl = .975*[10;.15;5;.2;55.6;3.87e-3;1e-6;4.5e-4;7.45e-4;5.22e-4;3e-6;3.3e-4;6e-9;5.37e-1;2.85e-1;7.79e-6;1e-6;4e-5;.01;.28;.05;.005;.005;.015;2.39;3e-4;.97];
xu = 1.025*[10;.15;5;.2;55.6;3.87e-3;1e-6;4.5e-4;7.45e-4;5.22e-4;3e-6;3.3e-4;6e-9;5.37e-1;2.85e-1;7.79e-6;1e-6;4e-5;.01;.28;.05;.005;.005;.015;2.39;3e-4;.97];

for ii = 1:N

    %Randomly sample parameters within acceptable ranges
    x = 2*rand(27,1) - 1;
    params = 1/2*(diag(xu - xl)*x + (xu + xl));
    A(ii,:) = x';
    
    %Create function handles
    f1 = @(X) params(1) + params(4)/(X(7) + params(5))*X(7)*X(1) - params(19)*X(1) - (params(6)*X(7) + params(7)*X(5))*X(1);
    f2 = @(X) params(27)*(params(6)*X(7) + params(7)*X(5))*X(1) + params(26)*X(3) - params(20)*X(2) - params(8)*X(2)*X(6);
    f3 = @(X) (1 - params(27))*(params(6)*X(7) + params(7)*X(5))*X(1) - params(26)*X(3) - params(21)*X(3);
    f4 = @(X) params(2) + params(9)*X(7)*X(4) - params(10)*X(7)*X(4) - params(22)*X(4);
    f5 = @(X) params(10)*X(7)*X(4) - params(23)*X(5) - params(11)*X(5)*X(6);
    f6 = @(X) params(3) + (params(12)*X(2) + params(13)*X(5))*X(6) - params(24)*X(6);
    f7 = @(X) params(14)*X(2) + params(15)*X(5) - params(16)*X(7)*X(1) - (params(17) + params(18))*X(7)*X(4) - params(25)*X(7);
    f = @(t,X) [f1(X),f2(X),f3(X),f4(X),f5(X),f6(X),f7(X)]';
    
    %Input initial conditions
    IC = [1000,0,0,30,0,500,1e-3];
    
    %Numerically solve system of ODEs
    options = odeset('events', @totalZero);
    [tout,fout] = ode23s(f,[0,2000],IC,options);
    q(ii) = fout(end,1);

end

%Calculate weights
A = [ones(N,1), A];
w = A\q;
w = w(2:end);
w = 1/norm(w).*w;

%Create plot
fig = figure(1);
plot(1:27,w,'.-b','MarkerSize',30);
title(['Algorithm 1.3 with ' int2str(N) ' Samples'],'Interpreter','latex','Fontsize',18,'FontWeight','bold')
xlabel('Parameters','Interpreter','latex','Fontsize',14)
ylabel('Parameter Weights','Interpreter','latex','Fontsize',14)
hold on
plot(linspace(0,30,100),.3*ones(100,1),'--k')
hold on
plot(linspace(0,30,100),-.3*ones(100,1),'--k')
ylim([-.75,.75])
xlim([0,28])
set(fig,'PaperUnits','inches','PaperSize',[10 8])
hgexport(fig, 'Alg1point3.pdf', hgexport('factorystyle'), 'Format', 'pdf');
