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

Question

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

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

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

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

我希望保留方程式插槽的这个不错的功能来切换 很容易在不同类型的功能反应之间。

这是可重现的示例：

1。定义模型

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。定义 ODE 系统

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。运行它

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

4。很好地绘制它

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。实施“稳态检查”

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

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。此外，您调用的是该函数的 external 版本，而不是对象，因此它当然无法访问“方程式”。

（希望是）好消息是 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)