# LOAD THE PACKAGES:
library(deSolve)
library(reshape2)
library(ggplot2)

# ENTRADAS DE MODELO:

# Vector que almacena el número inicial de personas en cada compartimento (en el paso 0)
#Toda la población que estamos modelando es susceptible a la infección.

initial_state_values <- c(S = 999,  # # ENTRADAS DE MODELO:

# Vector que almacena el número inicial de personas en cada compartimento (en el paso 0)
#Toda la población que estamos modelando es susceptible a la infección.
                          I = 1,       # La epidemia empieza con una persona infectada 
                          L = 0)       # no hay inmunidad previa en la población

## Vector que almacena los parámetros que describen las tasas de transición en unidades de días ^ -1
parameters <- c(betha= 0.002 ,  mu=0.05, ni= 0.02, rho=0.1 ,deltha= 0.001, alpha=  0.01 )   # La tasa de recuperación, que actúa sobre las personas infectadas.

# TIMESTEPS:

Vector almacenando la secuencia de pasos de tiempo para resolver el modelo en

de 0 a 60 días en intervalos diarios

times <- seq(from = 0, to = 300, by = 1)

FUNCIÓN MODELO SIR:

La función del modelo toma como argumentos de entrada (en el siguiente orden): tiempo, estado y parámetros

Frog_model <- function(time, state, parameters) {  

    with(as.list(c(state, parameters)), {  # dígale a R que desempaquete los nombres de variables de las entradas de estado y parámetros
        
     # Las ecuaciones diferenci ales
      dL <- betha * S -(1-L)*L - mu * L              
      # las personas salen del (-) compartimento S a una velocidad lambda (fuerza de infección)
      dS <- (1-L)*L - ni*S+ rho*I - deltha*S*I    # las personas se mueven hacia (+) las formas del compartimento I a una velocidad lambda, # y se mueven fuera de (-) el compartimento I a una velocidad gamma (recuperación)
      dI <- deltha*S*I-(ni+alpha+ rho)*I             # personas se mueven hacia (+) el compartimento R desde I a una velocidad gamma
      
     # Devuelve el número de personas en los compartimentos S, I y R en cada paso de tiempo
     # (en el mismo orden que las variables de estado de entrada)
    return(list(c(dL, dS, dI))) 
    })
  
}

#MODELO DE SALIDA (resolviendo las ecuaciones diferenciales):

Resolviendo las ecuaciones diferenciales usando el algoritmo de integración de odas

output_frog <- as.data.frame(ode(y = initial_state_values, 
                            times = times, 
                            func = Frog_model,
                            parms = parameters))
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 1.70952, R2 = 9.13807e-17
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 1.70952, R2 = 9.13807e-17
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 1.70952, R2 = 9.13807e-17
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 1.70952, R2 = 7.47845e-17
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 1.70952, R2 = 7.47845e-17
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 1.70952, R2 = 7.47845e-17
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 1.70952, R2 = 6.12025e-17
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 1.70952, R2 = 6.12025e-17
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 1.70952, R2 = 6.12025e-17
##  
## DLSODA-  Warning..Internal T (=R1) and H (=R2) are
##       such that in the machine, T + H = T on the next step  
##      (H = step size). Solver will continue anyway.
## In above message, R1 = 1.70952, R2 = 5.00871e-17
##  
## DLSODA-  Above warning has been issued I1 times.  
##      It will not be issued again for this problem.
## In above message, I1 = 10
##  
## DLSODA-  At current T (=R1), MXSTEP (=I1) steps   
##       taken on this call before reaching TOUT     
## In above message, I1 = 5000
##  
## In above message, R1 = 1.70952
## 
## Warning in lsoda(y, times, func, parms, ...): an excessive amount of work (>
## maxsteps ) was done, but integration was not successful - increase maxsteps
## Warning in lsoda(y, times, func, parms, ...): Returning early. Results are
## accurate, as far as they go
output_long <- melt(as.data.frame(output_frog), id = "time")               # convierte el conjunto de datos de salida en formato largo

ggplot(data = output_long,                                               # especificar objeto que contiene datos para trazar
       aes(x = time, y = value, colour = variable, group = variable)) +  # asignar columnas a ejes y grupos
  geom_line() +                                                          # representa los datos como líneas
  xlab("Time (days)")+                                                    # agregar etiqueta para el eje x
  ylab("Number of people") +                                              # agregar etiqueta para el eje y
  labs(colour = "Compartment")