%
  clear all

  global N

%
%  Set the initial conditions and constants. Load the data from SIR.txt.
%

  load SIR.txt

  S0 = 900;
  I0 = 100;
  R0 = 0;
  N = 1000;

  tf = 5;
  dt = 0.1;
  %t_data = 0:dt:tf;  % This is used when constructing the data
  t_data = SIR(:,1);
  I_data = SIR(:,2);
  n = length(t_data);
  p = 3;
  Y0 = [S0; I0; R0];
  delta = 0.2;
  gamma = 0.01;
  r = 0.8;
  h = 1e-8;
  q = [gamma,r,delta];
  
  sigma = 21.5;
  var = sigma^2;
  error = sigma*randn(size(t_data'));

%
% Numerically solve the system, construct data, and estimate the parameters.
%

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

  [t,Y] = ode45(@SIR_rhs,t_data,Y0,ode_options,q);
  I = Y(:,2);
%  I_data = I + error; % This is the data

  q_init = [0.05,0.5,0.5];
  modelfun = @(q)SIR_fun(q,Y0,t_data,I_data);
  [q_opt,fval] = fminsearch(modelfun,q_init);
  q_opt
  q

%
% Numerically solve the system with the perturbed parameters to compute the sensitivity matrix.
%

  for i=1:3
    if i==1
      q = [gamma+h,r,delta];
      [t,Ypert] = ode45(@SIR_rhs,t_data,Y0,ode_options,q);
      Ipert = Ypert(:,2);
      Sens(:,i) = (I - Ipert)/h; 
    elseif i==2
      q = [gamma,r+h,delta];
      [t,Ypert] = ode45(@SIR_rhs,t_data,Y0,ode_options,q);
      Ipert = Ypert(:,2);
      Sens(:,i) = (I - Ipert)/h;
    else 
      q = [gamma,r,delta+h];
      [t,Ypert] = ode45(@SIR_rhs,t_data,Y0,ode_options,q);
      Ipert = Ypert(:,2);
      Sens(:,i) = (I - Ipert)/h;
    end
  end
  
  Res = I - I_data;
  sigma2_est = (1/(n-p))*Res'*Res
  var 
  V = sigma2_est*inv(Sens'*Sens);
  rank(V)


%
%  Save the data.
%

  SIR_data = [t I_data];
%  save('SIR.txt','SIR_data','-ascii')

%
% Plot the results.
%

  

  figure(1)
  plot(t,Y(:,1),'k-',t,Y(:,2),'r--',t,Y(:,3),'b-.','linewidth',2)
  set(gca,'Fontsize',[20]);
  xlabel('Time (s)')
  ylabel('Responses S(t), I(t), R(t)')
  legend('Susceptible S(t)','Infectious I(t)','Recovered R(t)','Location','East')

  figure(2)
  plot(t,I,t,I_data,'x','linewidth',2)
  set(gca,'Fontsize',[20]);
  xlabel('Time (s)')
  ylabel('Infectious I(t)')
  legend('Model','Data','Location','NorthEast') 

  figure(3)
  plot(t,0*t,'k',t,Res,'x','linewidth',3)
  set(gca,'Fontsize',[20]);
  xlabel('Time (s)')
  ylabel('Residuals')