% Goal function for the optimization of beta and gamma 
% parameters of the SIR model for the Coronavirus epidemic 
% in Italy. 
% 
% To optimize, run in Matlab:
% opt = fmincon(@objective_function,[1.2 .95],[],[],[],[],[0 0],[])
%
% Caner Kazanci 4/15/20 (first version) 4/22/20 (updated) 

function err = objective_function( p )
format compact
data_I = [ 0.000000033333333
           0.000000033333333
           0.000000033333333
           0.000000033333333
           0.000000033333333
           0.000000033333333
           0.000000033333333
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000050000000
           0.000000333333333
           0.000001033333333
           0.000002583333333
           0.000003816666667
           0.000005366666667
           0.000007550000000
           0.000010916666667
           0.000014800000000
           0.000018800000000
           0.000028233333333
           0.000033933333333
           0.000041700000000
           0.000051483333333
           0.000064300000000
           0.000077266666667
           0.000098050000000
           0.000122916666667
           0.000152866666667
           0.000169150000000
           0.000207700000000
           0.000207700000000
           0.000294333333333
           0.000352616666667
           0.000412450000000
           0.000466333333333
           0.000525100000000
           0.000595216666667
           0.000683916666667
           0.000783683333333
           0.000892966666667
           0.000985633333333
           0.001065450000000
           0.001152933333333
           0.001239766666667
           0.001343150000000
           0.001441633333333
           0.001541200000000
           0.001628150000000
           0.001695650000000
           0.001763200000000
           0.001842900000000
           0.001920700000000
           0.001997116666667
           0.002077200000000
           0.002149133333333
           0.002209116666667
           0.002259766666667
           0.002323700000000
           0.002393766666667
           0.002459616666667
           0.002537850000000
           0.002606050000000
           0.002658600000000
           0.002708133333333
           0.002752583333333 ]';

data_t = 0:75;

time_range = 0:.01:75;
init = [1 3.3333e-08];
ode_func=@(t,x) [ -p(1)*x(1)*x(2); x(2)*(p(1)*x(1)-p(2)) ];
[t,x] = ode45( ode_func, time_range, init);

err = 0; 
for i = 1:76
  ind = find( abs(t - data_t(i)) < 1e-4 );
  err = err + ( x(ind,2) - data_I(i) )^2;
end
  
err = err * 1e6;

plot(t,x(:,2),'-'); hold on; 
plot(data_t,data_I,'*'); 
grid; 
xlabel('time'); 
ylabel('%infected'); hold off;
% axis([0 40 0 .00021])

fprintf('Error = %3.10f - Parameters = %2.1f, %2.1f\n',err,p)
