stopifnot(require(deSolve))

# Find the time derivatives of the Lotka-Volterra equations
  # Format fixed by the function ode() in the deSolve package
# Inputs: time (t) --- not used but required by the function ode()
  # two-dimensional state vector (state)
  # four-dimensional parameter vector (parameters)
# Output: a list containing the rates of change for the two state variables
# Presumes: state is a vector of length 2
  # Order of state is (lynx, hare)
  # parameters is a vector of length 4
  # Order of parameters is as in equations not included here
LotkaVolterra <- function(t, state, parameters) {
    dlynx <- (parameters[1]*state[2] - parameters[2])*state[1]
    dhare <- (parameters[3] - parameters[4]*state[1])*state[2]
    return(list(c(dlynx, dhare)))
}


# Solve the Lotka-Volterra equations from the fixed initial condition (1,1)
# with given parameters
# Inputs: four-dimensional parameters (lv.thetas)
  # times for outputs (times)
# Output: matrix with columns for times, relative number of lynxes, rleative number of hares
lv.latent <- function(lv.thetas, times=time(lynx)) {
    stopifnot(require(deSolve))
    initial.state <- c(1, 1)
    names(initial.state) <- c("lynx", "hare")
    return(ode(y=initial.state,
               times=times, func=LotkaVolterra,
               parms=c(lv.thetas)))
}

# Simulate observing only the predator numbers in a Lotka-Volterra predator-prey system, with noise
# Inputs: state-sequence data frame as returned by lv.latent (state.seq)
  # Observational-noise parameter (obs.theta)
# Output: vector of observed predator numbers
lv.observations <- function(state.seq, obs.theta) {
    return(rpois(n=nrow(state.seq), lambda=obs.theta*state.seq[,"lynx"]))
}


# Simulate the Lotka-Volterra model and return a time series of observed predator numbers
# Inputs: parameter vector (theta.vec)
  # time points for the time series (times)
# Output: a time series
lv.sim <- function(theta.vec, times=time(lynx)) {
    lv.out <- lv.latent(lv.thetas=theta.vec[1:4], times=times)
    lynx.series <- lv.observations(lv.out, obs.theta=theta.vec[5])
    return(ts(data=lynx.series, start=min(times), end=max(times)))
}

theta.ref <- c(0.94, 1.5, 0.33, 0.12, 550)