R simecol和稳态：rootfun无法访问SimObj方程槽

Question

我使用simecol包框架模拟生态模型（作为ODE系统），通过对象类SimObj（See here）非常容易地使用和共享生态模型。

我想实现一个稳态，一旦导数变得非常低，就会停止模拟。

根据这个vignette和这个example，你可以轻松实现它。

您只需要提供一个自定义求解器来检查导数的值。问题是自定义解算器看起来无法到达SimObj的equations插槽。

我想保持方程槽的这个很好的功能，以便在不同类型的功能响应之间轻松切换。

这是可重复的例子：

1. Define the model

library(simecol)

upca_model <- function() {
  new("odeModel",
    main = upca_ode,
    equations = list(
      f1 = function(x, y, k){x * y}, # Lotka-Volterra
      f2 = function(x, y, k){f1(x, y, k) / (1 + k * x)} # Holling II
      ),
    times = c(from = 0, to = 300, by = 0.1),
    parms = c(a = 1, b = 1, c = 10, alpha1 = 0.2, alpha2 = 1,
      k1 = 0.05, k2 = 0, wstar = 0.1),
    init = c(u = 10, v = 5, w = 0.1),
    solver = "lsoda"
    )
}

2. Define the ODE system

upca_ode <- function(time, init, parms) {
  u <- init["u"]
  v <- init["v"]
  w <- init["w"]
  with(as.list(parms), {
    du <- a * u - alpha1 * f(u, v, k1)
    dv <- -b * v + alpha1 * f(u, v, k1) - alpha2 * f(v, w, k2)
    dw <- -c * (w - wstar) + alpha2 * f(v, w, k2)
    list(c(du, dv, dw))
})
}

3. Run it

upca <- upca_model()
equations(upca)$f <- equations(upca)$f2
test <- sim(upca)

4. Plot it nicely

plotupca <- function(obj, ...) {
  o <- out(obj)
  matplot(o[, 1], o[, -1], type = "l", ...)
  legend("topright", legend = c("u", "v", "w"), lty = 1:3,, bg = "white",
    col = 1:3)
}

plotupca(test)

我们可以改变f方程，因此我们可以很容易地改变功能响应类型。

equations(upca)$f <- equations(upca)$f1
test <- sim(upca)
plotupca(test)

我们看到我们不需要长时间运行模拟，因为它似乎在大约100次步后达到稳定状态。

5. Implement the "steady-state checking"

因此，我们实现了一个解算器，一旦达到稳定状态就会停止模拟：

steady_state_upca <- function(time, init, func, parms) {
  root <- function(time, init, parms) {
    dstate <- unlist(upca_ode(time, init, parms))
    return(sum(abs(dstate)) - 1e-4)
  }
  lsodar(time, init, func, parms, rootfun = root)
}

equations(upca)$f <- equations(upca)$f1
solver(upca) <- steady_state_upca
test <- sim(upca)
#> Error in f(u, v, k1) : impossible de trouver la fonction "f"                

因此，找不到equation中定义的函数。

但是如果我在ODE系统中添加它，它就可以工作。

upca_ode <- function(time, init, parms) {
  u <- init["u"]
  v <- init["v"]
  w <- init["w"]
  # Definition of the function f:
  f <- function(x, y, k){x * y}
  with(as.list(parms), {
    du <- a * u - alpha1 * f(u, v, k1)
    dv <- -b * v + alpha1 * f(u, v, k1) - alpha2 * f(v, w, k2)
    dw <- -c * (w - wstar) + alpha2 * f(v, w, k2)

    list(c(du, dv, dw))
})
}

upca <- upca_model()
equations(upca)$f <- equations(upca)$f1
solver(upca) <- steady_state_upca
test <- sim(upca)
plotupca(test)

我们看到模拟已经提前停止（100而不是300），它已经停止，因为达到了稳定状态。

我的问题是：如何访问自定义求解器lsodar可访问的方程槽？

Answer 1

看了你的例子，我看到你自己的解算器有些“有线”，即不必要的复杂。解算器由sim函数调用并获得它所需的内容，因此无需在其中调用upca_ode。此外，您调用该函数的外部版本，而不是对象，因此它当然无法访问“方程式”。

（希望）这个好消息是包rootSolve包含已经可以使用的稳态求解器。

见下面的例子。下面的例子还表明，定义一个自己的求解器非常容易。请使用正确的参数顺序。

希望能帮助到你！

托马斯

library(simecol)
library(rootSolve)

upca <- new("odeModel",
            main = function(time, init, parms) {
              u      <- init[1]
              v      <- init[2]
              w      <- init[3]
              with(as.list(parms), {
                du <-  a * u           - alpha1 * f(u, v, k1)
                dv <- -b * v           + alpha1 * f(u, v, k1) +
                  - alpha2 * f(v, w, k2)
                dw <- -c * (w - wstar) + alpha2 * f(v, w, k2)
                list(c(du, dv, dw))
              })
            },
            equations  = list(
              f1 = function(x, y, k){x*y},           # Lotka-Volterra
              f2 = function(x, y, k){x*y / (1+k*x)}  # Holling II
            ),
            times  = c(from=0, to=300, by=0.1),
            parms  = c(a=1, b=1, c=10, alpha1=0.2, alpha2=1,
                       k1=0.05, k2=0, wstar=0.1),
            init = c(u=10, v=5, w=0.1),
            solver = "lsoda"
)
upca@equations$f <- upca@equations$f2

plotupca <- function(obj, ...) {
  o <- out(obj)
  matplot(o[, 1], o[, -1], type = "l", ...)
  legend("topright", legend = c("u", "v", "w"), lty = 1:3, bg = "white", col = 1:3)
}

test <- sim(upca)
plotupca(test)

upca@equations$f <- upca@equations$f1
test <- sim(upca)
plotupca(test)


steady_state <- function(y, times, func, parms, steady.method) {
  time <- if (steady.method == "stode") {
    0
  } else {
    c(times[1], Inf)
  }
  steady(y, time, func, parms, method=steady.method)$y
}

equations(upca)$f <- equations(upca)$f1
solver(upca) <- steady_state

# fast direct approach, does not work with all models
test <- sim(upca, steady.method="stode")
out(test)

## slower, simulates model until approx. steady
test <- sim(upca, steady.method="runsteady")
out(test)