stanode: R wrapper for modeling ODEs in Stan

#' stan_ode_generate
#' @description Generates a syntactically valid Stan program using the ODE system generated by
#' \code{\link{stan_lines}}.
#' @param obj Output of \code{stan_lines}.
#' @param path String path generated by \code{stan_prog}.
#' @param n_states Number of states in the ODE system.
#' @return A string with the Stan model which can be used in \code{\link[rstan]{stan}}. A more
#' readable representation of the model can be constructed by applying the
#' \code{\link[base]{cat}} function to the output.
#' @seealso \code{\link{stan_lines}}, \code{\link{stan_prog}}.
#' @examples
#' \dontrun{
#' # Define the ODE system as an R function
#' two_cpt <- function(t, y, parms) {
#'   with(as.list(c(y, parms)), {
#'     dy_gut = -ka * y_gut
#'     dy_cent = ka * y_gut - (CL/V_cent + Q/V_cent) * y_cent + (Q/V_peri) * y_peri
#'     dy_peri = (Q/V_cent) * y_cent - (Q/V_peri) * y_peri
#'     return(list(c(dy_gut=dy_gut, dy_cent=dy_cent, dy_peri=dy_peri)))
#'   })
#' }
#' # Create the time steps, inital state values, and parameter values
#' pars <- c("CL" = 10, Q = 13, "V_cent" = 20, "V_peri" = 73, ka = 3)
#' yini <- c("y_gut" = 0, "y_cent" = 0, "y_peri" = 0)
#' time_steps <- seq(0, 150, by = 0.01)
#' # Convert the R defined ODE into Stan syntax
#' sl <- stan_lines(two_cpt, state = yini,
#'                  pars = pars,
#'                  times = time_steps)
#' # Generate the useable Stan file
#' cat(stan_ode_generate(sl, has_events = TRUE, integrator = "rk45", n_states = 3))
#' }
#' @export

stan_ode_generate <- function(obj, path, n_states) {
  n_eqn <- length(obj)

  stan_wrapper <- readLines(path)

  sel <- which(stan_wrapper == "    #include \"user_func.stan\"")
  upr <- 1:(sel-1)
  lwr <- (sel+1):length(stan_wrapper)

  ode_eqns <- c("    // AUTO GENERATED CODE",
                paste0("    real dydt[",n_states,"];"),
                obj,
                "    return dydt;")
  out <- c(stan_wrapper[upr], ode_eqns, stan_wrapper[lwr])
  return(out)
}