#MIS 30.03.2020 zadatak 5 SIR model
# Zemlja: Japan

PopSt <- 126800000 # populacija japana

library(simecol)
library(ggplot2)

sir_model <- function(time, # current simulation time unit (day)
                      init, # number of susceptible, infected and recovered at this time unit
                      parms # beta and gamma values
) {
  # Function with unpacks given values for easier access
  # Note: with only works on lists, not on vectors
  with(as.list(c(init, parms)), {
    # Once unpacked we can now use these values by their names
    dS <- -beta * S * I
    dI <-  beta * S * I - gamma * I
    dR <-  gamma * I

    # This return is a result for with function
    return(list(c(dS, dI, dR)))
  })

  # Result from "with()" function is then returned by the main (sir_model(...)) function as we don't have to explicitly write return command when the last line of a function resolves in a value
}

param_values <- c(
  beta  = 0.000000355 / (PopSt / 1000000), # infection rate
  gamma = 0.2    # recovery rate
)

init_values <- c(
  S = PopSt-1,  # number of susceptibles at time = 0
  I =   1,  # number of infectious at time = 0
  R =   0   # number of recovered (and immune) at time = 0
)

times <- c(from=0, to=200, by=1)

SIR <- odeModel(main = sir_model,
                parms = param_values,
                times = times,
                init = init_values)

SIR <- sim(SIR)

SIR.data <- out(SIR)
head(SIR.data)

#graf
ggplot(SIR.data, aes(time)) +
  geom_line(aes(y = S, colour = "Podlozno zarazi")) +
  geom_line(aes(y = I, colour = "Zarazeno")) +
  geom_line(aes(y = R, colour = "Ozdravljeno")) +
  labs(color = "State")

R_0.func <- function(init, parms) {

  with(as.list(c(init, parms)), {
    # Total population * beta / gamma
    sum(S + I + R) * beta / gamma
  })

}

SIM.R_O <- R_0.func(init = init_values,
                    parms = param_values)

print(SIM.R_O)# R0 = 1.775

# -- { Prema sluzbenim podacima } -- #
# Ukupno slucajeva: 1,866
# Ukupno umrlih: 54
# Ukupno ozdravljeno: 424
# Prvi slucaj: Jan 14
# Proteklo dana od prvog slucaja: 76 dana ne uključujući današnji dan (30.3)

times2 <- c(from=0, to=76, by=1)

SIR <- odeModel(main = sir_model,
                parms = param_values,
                times = times2,
                init = init_values)

SIR <- sim(SIR)
SIR.data <- out(SIR)
print(SIR.data)

ggplot(SIR.data, aes(time)) +
  geom_line(aes(y = S, colour = "Podlozno zarazi")) +
  geom_line(aes(y = I, colour = "Zarazeno")) +
  geom_line(aes(y = R, colour = "Ozdravljeno")) +
  labs(color = "State")