Go back to Flu Mortality

library(deSolve)

dyn.load(“SIR_solver.so”)

data1 = c(106, 121, 136, 156, 183, 248, 279, 276, 239, 185, 160, 187, 150, 193, 153, 153, 138, 143, 119)

t1 = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18)

model_SIRn = function(t, theta) {

SIR = function(Time, State, Pars) {
with(as.list(c(State, Pars)), {

	DecreaseS = -1 * beta * S * I / N
	DecreaseI = -1 * alpha * I

	dS = DecreaseS
	dI = -1 * DecreaseS + DecreaseI

	return(list(c(dS, dI)))
	})
}	

with(as.list(theta), {
	
state1 = c(S = N - I1, I = I1)	
pars1 = c(beta = beta1, alpha = alpha1, N = N)
	
out1 = ode(y = state1, times = t, func = SIR, parms = pars1, 	method = 'ode45')
	
	return(-diff(out1[,"S"]))
})

}

model_SIRc = function(t, theta) {

SIR_solver = function(pars, t, state)
{
	ansS = double(1)
	
	   out = .C("SIR_solver",
                 params = as.double(pars),
                 t = as.double(t),
                 yin = as.double(state),
                 ansS = as.double(ansS))
                 
                 return(out$ansS)
}

with(as.list(theta), {
	
state1 = c(S = N - I1, I = I1)	
pars1 = c(beta = beta1, alpha = alpha1, N = N)

out1 = SIR_solver(pars1, t, state1)

	return(-diff(out1))
})

}

theta = c(beta1 = 0.3, alpha1 = 0.2, N = (sum(data1)*1.166), I1 = 1)

norm = model_SIRn(c(0:19), theta)

print(norm)

c = model_SIRc(c(0,19), theta)

print©