%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 for each time step
time = 1550:100:3350;  %Time values
h = 1e-6;   %Finite difference step size
trial = 1;  %Trial number (used when saving figures)

%Pre-allocate memory
q = zeros(N,numel(time));        %Output of interest (T-cell count after time(ii) days)
qplus = zeros(27,1);             %Perturbed output of interest
gradq = zeros(27,N,numel(time)); %Gradient of output of interest 
Xs = zeros(N,27,numel(time));    %To save the normalized paramters
w = zeros(27,numel(time));       %Weight vectors 
evalues = zeros(27,numel(time)); %Eigenvalues of the C matrix
diff = zeros(numel(time),1);     %Differences in largest and smallest element of gradq
I = eye(27);                     %27x27 identity matrix  

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

%Set initial conditions
IC = [1000,0,0,30,0,500,1e-3];

%Run simulation
for ii = 1:numel(time)
    ii
    for jj = 1:N
        
        %Randomly sample parameters within acceptable ranges
        Xs(jj,:,ii) = 2*rand(1,27) - 1;
        params = 1/2*(diag(xu - xl)*Xs(jj,:,ii)' + (xu + xl));
        
        %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)]';
        
        %Numerically solve system of ODEs
        options = odeset('events', @totalZero);
        [tout,fout] = ode23s(f,[0,time(ii)],IC,options);
        q(jj,ii) = fout(end,1);
        
        for kk = 1:27
            
            %Numerically solve perturbed system of ODEs
            xplus = Xs(jj,:,ii)' + h*I(:,kk);
            paramsplus = 1/2*(diag(xu - xl)*xplus + (xu + xl));
            f1 = @(X) paramsplus(1) + paramsplus(4)/(X(7) + paramsplus(5))*X(7)*X(1) - paramsplus(19)*X(1) - (paramsplus(6)*X(7) + paramsplus(7)*X(5))*X(1);
            f2 = @(X) paramsplus(27)*(paramsplus(6)*X(7) + paramsplus(7)*X(5))*X(1) + paramsplus(26)*X(3) - paramsplus(20)*X(2) - paramsplus(8)*X(2)*X(6);
            f3 = @(X) (1 - paramsplus(27))*(paramsplus(6)*X(7) + paramsplus(7)*X(5))*X(1) - paramsplus(26)*X(3) - paramsplus(21)*X(3);
            f4 = @(X) paramsplus(2) + paramsplus(9)*X(7)*X(4) - paramsplus(10)*X(7)*X(4) - paramsplus(22)*X(4);
            f5 = @(X) paramsplus(10)*X(7)*X(4) - paramsplus(23)*X(5) - paramsplus(11)*X(5)*X(6);
            f6 = @(X) paramsplus(3) + (paramsplus(12)*X(2) + paramsplus(13)*X(5))*X(6) - paramsplus(24)*X(6);
            f7 = @(X) paramsplus(14)*X(2) + paramsplus(15)*X(5) - paramsplus(16)*X(7)*X(1) - (paramsplus(17) + paramsplus(18))*X(7)*X(4) - paramsplus(25)*X(7);
            f = @(t,X) [f1(X),f2(X),f3(X),f4(X),f5(X),f6(X),f7(X)]';
            options = odeset('events', @totalZero);
            [tout,fout] = ode23s(f,[0,time(ii)],IC,options);
            qplus(kk) = fout(end,1);
            
        end
        
        %Calculate the gradients using finite differences
        gradq(:,jj,ii) = (qplus - q(jj,ii))/h;
        
    end
end

%Compute the weights, eigenvalues, and plot results
close all
for nn = 1:numel(time)
    
    %Compute the singular value decomposition of C
    [U,S,V] = svd(1/sqrt(N)*gradq(:,:,nn));
    w(:,nn) = U(:,1);
    w2 = U(:,2);
    
    %Compute the eigenvalues of C
    evalues(:,nn) = diag(S.^2);
    
    %Plot the eigenvalues of C on a log plot
    fig = figure;
    semilogy(1:27,evalues(:,nn),'.-b','MarkerSize',30)
    title(['Eigenvalues of C After ' int2str(time(nn)) ' Days (N = ' int2str(N) ')'],'Interpreter','latex','Fontsize',16,'FontWeight','bold','Position',[12.5 180 0])
    xlim([0,28])
    set(get(gca,'Title'),'Units','Normalized','Position',[.45,1.04])
    set(fig,'PaperUnits','inches','PaperSize',[10 8])
    hgexport(fig, ['Evalues' int2str(trial) '_' int2str(time(nn)) '.pdf'], hgexport('factorystyle'), 'Format', 'pdf');
 
    %Plot the weight vector
    fig = figure;
    plot(1:27,w(:,nn),'.-b','MarkerSize',30)
    title(['Weight Vector After ' int2str(time(nn)) ' Days (N = ' int2str(N) ')'],'Interpreter','latex','Fontsize',16,'FontWeight','bold','Position',[12.5 1.05 0])
    xlabel('Parameters','Interpreter','latex','Fontsize',14)
    ylabel('Parameter Weights','Interpreter','latex','Fontsize',14)
    set(get(gca,'Title'),'Units','Normalized','Position',[.45,1.04])
    xlim([0,28])
    ylim([-1,1])
    set(fig,'PaperUnits','inches','PaperSize',[10 8])
    hgexport(fig, ['WV' int2str(trial) '_' int2str(time(nn)) '.pdf'], hgexport('factorystyle'), 'Format', 'pdf');

    %Plot the absolute value of weight vector components
    fig = figure;
    plot(1:27,abs(w(:,nn)),'.-b','Markersize',30)
    title(['Parameter Weight After ' int2str(time(nn)) ' Days (N = ' int2str(N) ')'],'Interpreter','latex','Fontsize',16,'Fontweight','bold','Position',[12.5 1.05 0])
    xlabel('Parameters','Interpreter','latex','Fontsize',14)
    ylabel('Magnitude of Weight','Interpreter','latex','Fontsize',14)
    xlim([0,28])
    ylim([0,1])
    set(fig,'PaperUnits','inches','PaperSize',[10 8])
    hgexport(fig, ['WVmag' int2str(trial) '_' int2str(time(nn)) '.pdf'], hgexport('factorystyle'), 'Format', 'pdf');
    
    %Sufficient summary plot
    fig = figure;
    plot(Xs(:,:,nn)*w(:,nn),q(:,nn),'ko');
    title(['Sufficient Summary Plot After ' int2str(time(nn)) ' Days (N = ' int2str(N) ')'],'Interpreter','latex','Fontsize',16,'FontWeight','bold','Position',[-.15 1130 0])
    xlabel('$w^T x_j$','Interpreter','latex','FontSize',14)
    ylabel('T-cell count','Interpreter','latex','FontSize',14)
    set(get(gca,'Title'),'Units','Normalized','Position',[.45,1.04])
    xlim([-2,2])
    ylim([0,1100])
    axis square;
    grid on;
    set(fig,'PaperUnits','inches','PaperSize',[10 8])
    hgexport(fig, ['SSP' int2str(trial) '_' int2str(time(nn)) '.pdf'], hgexport('factorystyle'), 'Format', 'pdf');

    %Find the difference of max and min gradq to check for errors
    diff(nn) = max(max(gradq(:,:,nn))) - min(min(gradq(:,:,nn)));
   
end

%Save the trial data
%save(['trial' int2str(trial) 'data.mat'],'evalues','Xs','q','w','diff','time')