%
%                        Example8_25_6param.m 
%

  clear all

  global coef Y0 wt

%
%  Set the number of samples N, parameters and initial conditions
%
%  Employ N = 1000 samples for final figures but N = 100 for debugging.
%
%  Note: We split the coefficients between params, which is passed to ode45
%       and coef which is global.  The former comprises the set employed in 
%       Chapters 12 and 13 for Bayesian inference and uncertainty propagation.
%

  N = 1000;

  lam1 = 1e+4;
  d1 = .01;
  epsilon = 0;
  k1 = 8.0e-7;
  lam2 = 31.98;
  d2 = 0.01;
  f = 0.34;
  k2 = 1e-4;
  delta = 0.7;
  m1 = 1.0e-5;
  m2 = 1.0e-5;
  NT = 100;
  c = 13;
  rho1 = 1;
  rho2 = 1;
  lamE = 1;
  bE = 0.3;
  Kb = 100;
  dE = 0.25;
  Kd = 500;
  deltaE = 0.1;

  coef = [epsilon,k1,lam2,d2,f,NT,c,rho1,rho2,lamE,dE,deltaE,m1,m2,Kd];
  parameters = [bE; delta; d1; k2; lam1; Kb];

  n = 41;
  p = 6;
  Y0 = [0.9e+6; 4000; 1e-1; 1e-1; 1; 12];

%
% Compute the time stpes and Effector cell population on the interval [0,200]. 
%

  tfinal = 200;
  dt = 0.1; 
  t_data = 0:dt:tfinal;
  N_data = length(t_data);
  ode_options = odeset('RelTol',1e-7);
  [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
  E = Y(:,6);

%
% Compute the local sensitivity matrix and Fisher information matrix for parameter
% subset selection.  We consider the parameters theta = [bE, delta, d1, k2, lam1, Kb].  
% We consider responses comprised of E(t), and [T_1(t), E(t)].
%

  h = 1e-12;

  coef = [epsilon,k1,lam2,d2,f,NT,c,rho1,rho2,lamE,dE,deltaE,m1,m2,Kd];

  bE_complex = complex(bE,h);
  parameters = [bE_complex; delta; d1; k2; lam1; Kb];
  [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
  %E_bE = imag(Y(:,6))/h;
  E_bE =  [imag(Y(:,1))/h; imag(Y(:,2))/h; imag(Y(:,3))/h; imag(Y(:,4))/h; imag(Y(:,5))/h; imag(Y(:,6))/h]; 

  delta_complex = complex(delta,h);
  parameters = [bE; delta_complex; d1; k2; lam1; Kb];
  [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
  %E_delta = imag(Y(:,6))/h;
  E_delta =  [imag(Y(:,1))/h; imag(Y(:,2))/h; imag(Y(:,3))/h; imag(Y(:,4))/h; imag(Y(:,5))/h; imag(Y(:,6))/h];

  d1_complex = complex(d1,h);
  parameters = [bE; delta; d1_complex; k2; lam1; Kb];
  [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
  %E_d1 = imag(Y(:,6))/h;
  E_d1 =  [imag(Y(:,1))/h; imag(Y(:,2))/h; imag(Y(:,3))/h; imag(Y(:,4))/h; imag(Y(:,5))/h; imag(Y(:,6))/h];

  k2_complex = complex(k2,h);
  parameters = [bE; delta; d1; k2_complex; lam1; Kb];
  [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
  %E_k2 = imag(Y(:,6))/h;
  E_k2 =  [imag(Y(:,1))/h; imag(Y(:,2))/h; imag(Y(:,3))/h; imag(Y(:,4))/h; imag(Y(:,5))/h; imag(Y(:,6))/h];

  lam1_complex = complex(lam1,h);
  parameters = [bE; delta; d1; k2; lam1_complex; Kb];
  [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
  %E_lam1 = imag(Y(:,6))/h;
  E_lam1 =  [imag(Y(:,1))/h; imag(Y(:,2))/h; imag(Y(:,3))/h; imag(Y(:,4))/h; imag(Y(:,5))/h; imag(Y(:,6))/h];

  Kb_complex = complex(Kb,h);
  parameters = [bE; delta; d1; k2; lam1; Kb_complex];
  [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
  %E_Kb = imag(Y(:,6))/h;
  E_Kb =  [imag(Y(:,1))/h; imag(Y(:,2))/h; imag(Y(:,3))/h; imag(Y(:,4))/h; imag(Y(:,5))/h; imag(Y(:,6))/h];


  Sens = [bE*E_bE delta*E_delta d1*E_d1 k2*E_k2 lam1*E_lam1 Kb*E_Kb];
  Fisher = Sens'*Sens;

  [V_f,D_f] = eig(Fisher);

  [~,Svals,V] = svd(Sens,'econ');
  Svals
  (Svals(6,6)/Svals(1,1))^2
  V(:,6) 

%
% Randomly sample parameters
%

  per = 0.1;
  per2 = 0.7;
  bE_param = bE*((1-per) + 2*per*rand(N,1));
  delta_param = delta*((1-per) + 2*per*rand(N,1));
  d1_param = d1*((1-per) + 2*per*rand(N,1));
  k2_param = k2*((1-per) + 2*per*rand(N,1));
  lam1_param = lam1*((1-per) + 2*per*rand(N,1));
  Kb_param = Kb*((1-per2) + 2*per2*rand(N,1));

%  bE_param = lognrnd(log(bE),0.05*bE,N,1);
  delta_param = lognrnd(log(delta),0.02*delta,N,1);
  d1_param = lognrnd(log(d1),0.02*d1,N,1);
  k2_param = lognrnd(log(k2),0.02*k2,N,1);
%  lam1_param = lognrnd(log(lam1),0.05*lam1,N,1);
  Kb_param = lognrnd(log(Kb),0.02*Kb,N,1);

  for j=1:N
    parameters = [bE_param(j); delta_param(j); d1_param(j); k2_param(j); lam1_param(j); Kb_param(j)];
    [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
    E1(j) = Y(N_data,6);
    E_T1(j) = Y(N_data,1);
  end

  for j=1:N
    parameters = [bE_param(j); delta_param(j); d1_param(j); k2_param(j); lam1_param(j); Kb];
    [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
    E2(j) = Y(N_data,6);
    E_T2(j) = Y(N_data,1);
  end

  xvals = 0:0.1:100;
  xvals2 = 1.62e+5:10:1.65e+5;
  dens1 = ksdensity(E1,xvals);
  dens2 = ksdensity(E2,xvals);
  dens3 = ksdensity(E_T1,xvals2);
  dens4 = ksdensity(E_T2,xvals2);

%
% Plot the distributions for E1 and E2
%

  figure(1)
  plot(xvals,dens1,'b',xvals,dens2,'b--','linewidth',3)
  axis([0 100 0 0.05])
  set(gca,'Fontsize',[24]);
  xlabel('E(200,\theta)')
  ylabel('PDF') 
  legend('All Random', 'Fixed K_b', 'Location','NorthEast')
  %print -depsc Fixed_Kb

  figure(2)
  plot(xvals2,dens3,'b-',xvals2,dens4,'b--','linewidth',3)
  axis([1.62e+5 1.65e+5 0 1.4e-3])
  set(gca,'Fontsize',[24]);
  xlabel('T_1(200,\theta')
  ylabel('PDF')
  legend('All Random', 'Fixed K_b', 'Location','NorthEast')


%
%  End hiv_pss
%