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

Question

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

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

据此 vignette and this 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 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

看了你的例子后，我发现你自己的求解器有点"wired"，即不必要的复杂。求解器由 sim 函数调用并得到它需要的东西，所以不需要在里面调用 upca_ode。此外，您调用该函数的 external 版本，而不是对象，因此它当然无法访问 "equations".

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

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

R simecol and steady-state: rootfun cannot access to the SimObj equations slot

simulation

functional-programming

r

ode

1。定义模型

2。定义 ODE 系统

3。运行它

4。画得好

5。实施 "steady-state checking"

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

R simecol and steady-state: rootfun cannot access to the SimObj equations slot

simulation

functional-programming

r

ode

1。定义模型

2。定义 ODE 系统

3。 运行它

4。画得好

5。实施 "steady-state checking"

3。运行它