%
%                      Example8_25_2param.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 Effector cell population on the interval [0,50].
%

  t_data = 0:0.1:50;
  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 = [rho1, lamE].
%

  h = 1e-12;

  rho1_complex = complex(rho1,h);
  coef = [epsilon,k1,lam2,d2,f,NT,c,rho1_complex,rho2,lamE,dE,deltaE,m1,m2,Kd];
  parameters = [bE; delta; d1; k2; lam1; Kb];
  [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
  E_rho1 = imag(Y(:,6))/h;

  lamE_complex = complex(lamE,h);
  coef = [epsilon,k1,lam2,d2,f,NT,c,rho1,rho2,lamE_complex,dE,deltaE,m1,m2,Kd];
  parameters = [bE; delta; d1; k2; lam1; Kb];
  [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
  E_lamE = imag(Y(:,6))/h;

  Sens = [E_rho1 E_lamE];
  Fisher = Sens'*Sens;

  [V_f,D_f] = eig(Fisher);

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

%
% Randomly sample parameters at N values specified at beginning of code.
%

  per = 0.1;
  rho1_param = rho1*((1-per) + 2*per*rand(N,1));
  lamE_param = lamE*((1-per) + 2*per*rand(N,1));

  for j=1:N
    coef = [epsilon,k1,lam2,d2,f,NT,c,rho1_param(j),rho2,lamE_param(j),dE,deltaE,m1,m2,Kd];
    parameters = [bE; delta; d1; k2; lam1; Kb];
    [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
    E1(j) = Y(N_data,6);
  end

  for j=1:N
    coef = [epsilon,k1,lam2,d2,f,NT,c,rho1,rho2,lamE_param(j),dE,deltaE,m1,m2,Kd];
    parameters = [bE; delta; d1; k2; lam1; Kb];
    [t,Y] = ode15s(@hiv_rhs,t_data,Y0,ode_options,parameters);
    E2(j) = Y(N_data,6);
  end

  xvals = 30:0.1:45; 
  dens1 = ksdensity(E1,xvals);
  dens2 = ksdensity(E2,xvals);
%
% Plot the distributions for E1 and E2
%

  figure(1)
  plot(xvals,dens1,'b',xvals,dens2,'b--','linewidth',3)
  axis([32 44 0 0.2])
  set(gca,'Fontsize',[24]);
  xlabel('E(50,\theta)')
  ylabel('PDF') 
  legend('Random \rho_1, \lambda_E', 'Random \lambda_E, Fixed \rho_1', 'Location','NorthEast')
%  print -depsc E_dist_2param

%
%  End Example8_25_2param.m
%