%options = statset('Display', 'iter', 'TolFun', 1e-19, 'MaxIter',200, 'TolX', 1e-19);
clear all
global y N Icc1 t1
N = 2000000;

t1=1:69;
Icc0 = [554	771	1208	1870	2613 4349	5739	7417	9308	11289	13748	16369	19383	22942	26302	28985	31774	33738	35982	37626	38791	51591	55748	56873	57416	57934	58016	57805	56301	54921	53284	52093	49824	47765	45600	43258	39919	37414	35129	32616	30004	27423	25353	23784	22179	20533	19016	17721	16136	14831	13524	12088	10733	9893	8967	8056	7263	6569	6013	5353	5120	4735	4287	3947	3460	3128	2691	2691	2161];
%N=sum(Icc0);
Icc1=Icc0/N;

guess =  [0.001 0.001 0.001 0.00 0.00 0.00];
A = [];
b = [];
Aeq = [];
beq = [];
lb = [0.01 0.001 0.001 0.001 0.001 0.001];
ub = [1.5 1.5 0.01 1.0 0.25 0.6];


[p,error] = fmincon(@DBfit, guess, A, b, Aeq, beq, lb, ub)

plot(t1,y(:,3)*N,t1,Icc1*N,'x', 'MarkerSize',6, 'LineWidth',2.50,'Color','black')
%title('','FontSize', 16)
xlabel('Time (days)','FontSize', 16)
ylabel('Active COVID-19 cases: China','FontSize', 16)
legend('Model estimate ',' Real data','Location','Best');

sqrterror = sqrt(error)


function dy = DAfit(t,y,p)

mu = 0.002; % Chen et al. Infectious Diseases of Poverty 
delta  = 0.113; %.High estimated per day
gamma  = 0.172; % Chen et al. Infectious Diseases of Poverty
xi = 0.100;   % Chen et al. Infectious Diseases of Poverty
 
dy = [mu-(1-p(6))*p(1)*y(1)*y(3)-p(2)*y(1)*y(6)-(mu+p(3))*y(1); 
(1-p(6))*p(1)*y(1)*y(3)+p(2)*y(1)*y(6)-(mu+p(4))*y(2); 
p(4)*y(2)-(mu+p(5)+delta+gamma)*y(3); 
p(5)*y(3)-(mu+delta+gamma)*y(4); 
p(3)*y(1)+(delta+gamma)*(y(3)+y(4))-mu*y(5); 
xi*(y(3)-y(6))];

end

function value = DBfit(p)

global y N Icc1 t1

s0 =0.9;
e0 = 0.0;
i0 = Icc1(1);
q0 = 0.0;
r0 = 0.1;
w0 = 0.01;

II = [s0 e0 i0 q0 r0 w0];

[t, y] = ode45(@DAfit, t1, II, [], p);
size(y);
error = (Icc1'-y(:,3)).^2;
value = sum(error);

end