%
%                    Example11_24.m
%
%  Requires: SIR_fun2.m, SIR_rhs2.m SIR2.txt, SIR2_integer.txt
%

  clear all

  global N

%
%  Load data which consists of 141 points of infectious data along with the
%  associated times.
%

  %load SIR2.txt
  load SIR2_integer.txt

%
%  Set parameters and initial conditions.  These values are close to the flu data
%  but with larger I0 and smaller variance to prevent negative values of I(t).
%

  S0 = 760;
  I0 = 12;
  R0 = 0;
  N = S0 + I0 + R0;

  beta = 0.0022;         % Estimated values for flu data: beta = 2.3040e-3;
  gamma = 0.4470;       %                                 gamma = 0.4284;
  beta_mean = beta;
  gamma_mean = gamma;
  params = [beta gamma];

  h = 1e-16;
  tf = 14;
  dt = 0.1;
  t_data = SIR2_integer(:,1);
%  I_data = SIR2(:,2);
  I_data = SIR2_integer(:,2);
  t_data = 0:dt:tf;
  Y0 = [S0; I0];

  n = length(t_data);
  p = 2;

%
% Specify the parameters and compute the solution I(t).  When integrating
% the system, we take advantage of the fact that R(t) = N - S(t) - I(t).
%
% Construct synthetic data having the observation variance sigma^2 = 81;
% Note that this is significantly less than the observation variance
% sigma^2 = 322.8 estimated for the flu data.
%

  q_init = [beta gamma];
  ode_options = odeset('RelTol',1e-8);
  [t,Y] = ode45(@SIR_rhs2,t_data,Y0,ode_options,q_init);
  Y(:,3) = N - Y(:,1) - Y(:,2);
  I = Y(:,2);

  sigma = 9;
  sigma2 = sigma^2;
  error = sigma*randn(size(t_data'));
%  I_data = I + error;
%  I_data_integer = round(I_data);
%  SIR_data = [t_data' I_data];
%  SIR_data_integer = [t_data' I_data_integer];
%  save('SIR2.txt','SIR_data','-ascii')
%  save('SIR2_integer.txt','SIR_data_integer','-ascii')

  Res = I_data - I;

%
%  Plot the data and model
%

  figure(1)
  plot(t,Y(:,1),'k-',t,Y(:,2),'r--',t,Y(:,3),'b-.','linewidth',3)
  axis([0 tf 0 800])
  set(gca,'Fontsize',[22]);
  xlabel('Time (Days)')
  ylabel('Responses S(t), I(t), R(t)')
  legend('Susceptible S(t)','Infectious I(t)','Recovered R(t)','Location','East')
  %print -depsc SIR2_responses

  figure(2)
  plot(t,I,'k-',t,I_data,'bx',t,0*t,'k-','linewidth',3)
  axis([0 tf 0 350])
  set(gca,'Fontsize',[22]);
  xlabel('Time (Days)')
  ylabel('Infectious I(t)')
  legend('Model','Data','Location','NorthEast')
  %print -depsc SIR2_model_data

  figure(3)
  plot(t,Res,'xb',t,0*t,'k-','linewidth',3)
  axis([0 tf -25 25]);
  hold on
  set(gca,'Fontsize',[22]);
  plot(t,2*sigma*ones(length(t)),'--k','linewidth',3)
  plot(t,-2*sigma*ones(length(t)),'--k','linewidth',3)
  hold off
  xlabel('Time (Days)')
  ylabel('Residuals')
  %print -depsc SIR2_residuals


%
% Estimate the parameters beta and gamma and compute the resulting trajectories I(t)
%

  ode_options = odeset('RelTol',1e-6);

  q_init = [beta gamma];
  modelfun = @(q)SIR_fun2(q,Y0,t_data,I_data);
  [q_opt,fval] = fminsearch(modelfun,q_init);

  beta_opt = q_opt(1)
  gamma_opt = q_opt(2)

  [t,Y] = ode45(@SIR_rhs2,t_data,Y0,ode_options,q_opt);
  Y(:,3) = N - Y(:,1) - Y(:,2);
  I = Y(:,2);

%
% Use the complex-step derivative approximation to compute the sensitivity matrix. Compute
% the observation error estimate and V.
%

  Y0 = [S0; I0];
  for i=1:p
    if i==1
      beta_complex = complex(beta_opt,h);
      q = [beta_complex,gamma];
      [t,Ypert] = ode45(@SIR_rhs2,t_data,Y0,ode_options,q);
      I_beta = imag(Ypert(:,2))/h;
      Sens(:,i) = I_beta;
    else
      gamma_complex = complex(gamma_opt,h);
      q = [beta,gamma_complex];
      [t,Ypert] = ode45(@SIR_rhs2,t_data,Y0,ode_options,q);
      I_gamma = imag(Ypert(:,2))/h;
      Sens(:,i) = I_gamma;
    end
  end
  
  Res = I_data - I;
  sigma2_est = (1/(n-p))*Res'*Res
  V = sigma2_est*inv(Sens'*Sens)

%%
% Construct the 95% confidence intervals.
%

  alpha = 0.05;
  
  tval = tinv(1-alpha/2,n-p); 
  int_beta = [beta_opt - sqrt(V(1,1))*tval  beta_opt + sqrt(V(1,1))*tval]
  int_r = [gamma_opt - sqrt(V(2,2))*tval  gamma_opt + sqrt(V(2,2))*tval]

%
% 
