%
%                      Example12_21.m 
%

  clear all

  global coef Y0 wt

%
%  Set parameters and initial conditions
%

  %load hiv_data_prev
  load hiv_data
  t_data = hiv_data(:,1);
  T1_data = hiv_data(:,2);
  T2_data = hiv_data(:,3);
  T1s_data = hiv_data(:,4);
  T2s_data = hiv_data(:,5);
  V_data = hiv_data(:,6);
  E_data = hiv_data(:,7);

  data.xdata = hiv_data(:,1);
  data.ydata = hiv_data(:,2:7);

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

  %coef = [lam1,epsilon,k1,lam2,d2,f,m1,m2,NT,c,rho1,rho2,lamE,Kb,dE,Kd,deltaE];
  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 residual and initial error variance.
%

  %wt = [1 1e+2 1 1e+2 1 1e+3];
  wt = [1 1 1 1 1 1];

  ode_options = odeset('RelTol',1e-6);
  [t,Y] = ode15s(@hiv_rhs_dram,t_data,Y0,ode_options,parameters);

  T1 = Y(:,1);
  T2 = Y(:,2);
  T1s = Y(:,3);
  T2s = Y(:,4);
  V = Y(:,5);
  E = Y(:,6);
  
  %sigma = .05;
  %error = sigma*Y.*randn(size(Y));
  %Y_data = Y + error;
  %T1_data = Y_data(:,1);
  %T2_data = Y_data(:,2);
  %T1s_data = Y_data(:,3);
  %T2s_data = Y_data(:,4);
  %V_data = Y_data(:,5);
  %E_data = Y_data(:,6);
  
  
  figure(1)
  semilogy(t_data,T1,t_data,T1_data,'--')
  
  figure(2)
  semilogy(t_data,T2,t_data,T2_data,'--')
  
  figure(3)
  semilogy(t_data,T1s,t_data,T1s_data,'--')
   
  figure(4)
  semilogy(t_data,T2s,t_data,T2s_data,'--')
  
  figure(5)
  semilogy(t_data,V,t_data,V_data,'--')
  
  figure(6)
  semilogy(t_data,E,t_data,E_data,'--')
  
  close all

  Res1 = wt(1)*(T1-T1_data);
  Res2 = wt(2)*(T2-T2_data); 
  Res3 = wt(3)*(T1s-T1s_data)
  Res4 = wt(4)*(T2s-T2s_data); 
  Res5 = wt(5)*(V-V_data);
  Res6 = wt(6)*(E-E_data);

  ss1 = Res1'*Res1;
  sigma1 = ss1/(n-p); 
  ss2 = Res2'*Res2;
  sigma2 = ss2/(n-p);
  ss3 = Res3'*Res3;
  sigma3 = ss3/(n-p);
  ss4 = Res4'*Res4;
  sigma4 = ss4/(n-p);
  ss5 = Res5'*Res5;
  sigma5 = ss5/(n-p);
  ss6 = Res6'*Res6;
  sigma6 = ss6/(n-p);
 
%
% Set up the model components for the DRAM implementation. 
%

  clear model options

  tmin = [bE; delta; d1; k2; lam1; Kb];

  params = {
    {'bE', tmin(1), 0}
    {'delta', tmin(2), 0}
    {'d1', tmin(3), 0}
    {'k2', tmin(4), 0}
    {'lam1', tmin(5), 0}
    {'Kb', tmin(6), 0}};
    
  model.ssfun = @hivss;
  model.S20 = [sigma1 sigma2 sigma3 sigma4 sigma5 sigma6];
  model.N0 = [1 1 1 1 1 1];
%  model.S20 = [sigma1 sigma5];
%  model.N0 = [1 1];
%  model.S20 = sigma6;
%  model.N0 = 1;
  options.nsimu = 30000;   % March 7   Change to 20000
%  options.nsimu = 3000;
  options.updatesigma = 1;

  [results,chain,s2chain] = mcmcrun(model,data,params,options);

  figure(1); clf
  mcmcplot(chain,[],results,'chainpanel');

  figure(2); clf
  mcmcplot(chain,[],results,'pairs');

  figure(3); clf
  mcmcplot_custom(chain,[],results,'denspanel',1); 

  figure(4); clf
  mcmcplot(chain,[],results,'hist',20);

%
% Process the chains to plot the densities.
%

  bE_vals = chain(:,1);
  delta_vals = chain(:,2);
  d1_vals = chain(:,3);
  k2_vals = chain(:,4);
  range_bE = max(bE_vals) - min(bE_vals);
  range_delta = max(delta_vals) - min(delta_vals);
  range_d1 = max(d1_vals) - min(d1_vals);
  range_k2 = max(k2_vals) - min(k2_vals);
  bE_min = min(bE_vals)-range_bE/10;
  bE_max = max(bE)+range_bE/10;
  delta_min = min(delta_vals)-range_delta/10;
  delta_max = max(delta_vals)+range_delta/10;
  d1_min = min(d1_vals)-range_d1/10;
  d1_max = max(d1)+range_d1/10;
  k2_min = min(k2_vals)-range_k2/10;
  k2_max = max(k2_vals)+range_k2/10;
  [bandwidth_bE,density_bE,bE_mesh,cdf_bE]=kde(bE_vals);
  [bandwidth_delta,density_delta,delta_mesh,cdf_delta]=kde(delta_vals);
  [bandwidth_d1,density_d1,d1_mesh,cdf_d1]=kde(d1_vals);
  [bandwidth_k2,density_k2,k2_mesh,cdf_k2]=kde(k2_vals);

  figure(11)
  plot(bE_mesh,density_bE,'k-','linewidth',3)
%  axis([-19.5 -17.5 0 3])
  set(gca,'Fontsize',[22]);
  xlabel('Parameter b_E')

  figure(12)
  plot(delta_mesh,density_delta,'k-','linewidth',3)
%  axis([1.8e-3 2e-3 0 3e4])
  set(gca,'Fontsize',[22]);
  xlabel('Parameter \delta')

  figure(13)
  plot(d1_mesh,density_d1,'k-','linewidth',3)
%  axis([-19.5 -17.5 0 3])
  set(gca,'Fontsize',[22]);
  xlabel('Parameter d_1')

  figure(14)
  plot(k2_mesh,density_k2,'k-','linewidth',3)
%  axis([1.8e-3 2e-3 0 3e4])
  set(gca,'Fontsize',[22]);
  xlabel('Parameter k_2')

%
% Construct the confidence limits for predictions.
%

  figure(20)
  t = linspace(0,200)';
  out = mcmcpred(results,chain,s2chain,t,@hiv_fun,5000)
  mcmcpredplot_custom(out)

%  subplot(6,1,1)
%  axis([0 200 0 2e+6])
%  subplot(6,1,2)
%  axis([0 200 0 5000])
%  subplot(6,1,3)
%  axis([0 200 0 5e+6])
%  subplot(6,1,2)
%  axis([0 200 0 5000])
%  subplot(6,1,4)
%  axis([0 200 0 2000])
%  subplot(6,1,5)
%  axis([0 200 0 2e+6])
%  subplot(6,1,6)
%  axis([0 200 0 50])

%  hold on
%  for i = 1:6
%     subplot(6,1,i)
%     hold
%     plot(t_data,data.ydata(:,i),'s')
%     hold off
%  end
  

%
%  End hiv_dram
%