##  https://strimas.com/post/lotka-volterra/

library(tidyverse)
library(deSolve)
library(FME)

params1 <- c(alpha=1.2,beta=0.2, delta =0.5, gamma=0.2)
params2 <- c(alpha=1.0,beta=0.2, delta =0.5, gamma=0.2)
params3 <- c(alpha=0.8,beta=0.2, delta =0.5, gamma=0.2)

init <- c(x=1, y=1)

times <- seq(0,100,by=1)


deriv <- function(t, state, pars) {
  with(as.list(c(state, pars)), {
    d_x <- alpha * x - beta * x * y
    d_y <- delta * beta * x * y - gamma * y
    return(list(c(x = d_x, y = d_y)))
  })
}
lv_results1 <- ode(init, times, deriv, params1)
lv_results2 <- ode(init, times, deriv, params2)
lv_results3 <- ode(init, times, deriv, params3)


lv_model <- function(pars, times = seq(0, 50, by = 1)) {
  # initial state 
  state <- c(x = 1, y = 2)
  # derivative
  deriv <- function(t, state, pars) {
    with(as.list(c(state, pars)), {
      d_x <- alpha * x - beta * x * y
      d_y <- delta * beta * x * y - gamma * y
      return(list(c(x = d_x, y = d_y)))
    })
  }
  # solve
  ode(y = state, times = times, func = deriv, parms = pars)
}

lv_results1 <- lv_model(pars = params1, times = seq(0, 50, by = 0.25)) 
lv_results2 <- lv_model(pars = params2, times = seq(0, 50, by = 0.25)) 
lv_results3 <- lv_model(pars = params3, times = seq(0, 50, by = 0.25)) 

lv_results1 %>% 
  data.frame() %>% 
  gather(var, pop, -time) %>% 
  mutate(var = if_else(var == "x", "Prey", "Predator")) %>% 
  ggplot(aes(x = time, y = pop)) +
  geom_line(aes(color = var)) +
  scale_color_brewer(NULL, palette = "Set1") +
  labs(title = "Lotka-Volterra predator prey model",
       subtitle = paste(names(params1), params1, sep = " = ", collapse = "; "),
       x = "Time", y = "Population density")  


plot(lv_results3[,2],lv_results3[,3],xlab="Prey",ylab="Predator",ylim=c(0,20),xlim=c(0,20),
     type="l",col="blue",cex.axis=1.5,cex.lab=1.5)
lines(lv_results2[,2],lv_results2[,3],col="red")
lines(lv_results1[,2],lv_results1[,3],col="green")