r - 编译的 C ODE 使用 deSolve 为 R 提供不同的结果

Question

我有一个 ODE，我想使用从 R 的 deSolve 包中调用的编译 C 代码来解决它。有问题的 ODE 是指数衰减模型 (y'=-d* exp(g* time)*y):但是从 R 中运行编译后的代码会给 R 的 native deSolve 带来不同的结果。就像它们被翻转 180º 一样。怎么回事？

C代码实现

/* file testODE.c */
#include <R.h>
static double parms[4];
#define C parms[0] /* left here on purpose */
#define d parms[1]
#define g parms[2]

/* initializer  */
void initmod(void (* odeparms)(int *, double *))
{
  int N=3;
  odeparms(&N, parms);
}

/* Derivatives and 1 output variable */
void derivs (int *neq, double t, double *y, double *ydot,
             double *yout, int *ip)
{
  // if (ip[0] <1) error("nout should be at least 1");
  ydot[0] = -d*exp(-g*t)*y[0];
}
/* END file testODEod.c */

R 实现 - Native deSolve

  testODE <- function(time_space, initial_contamination, parameters){
  with(
    as.list(c(initial_contamination, parameters)),{
      dContamination <- -d*exp(-g*time_space)*Contamination
      return(list(dContamination))
    }
  )
}

parameters <- c(C = -8/3, d = -10, g =  28)
Y=c(y=1200)
times <- seq(0, 6, by = 0.01)
initial_contamination=c(Contamination=1200) 
out <- ode(initial_contamination, times, testODE, parameters, method = "radau",atol = 1e-4, rtol = 1e-4)

plot(out)

R 实现 - 从 deSolve 运行编译代码

library(deSolve)
library(scatterplot3d)
dyn.load("Code/testODE.so")

Y <-c(y1=initial_contamination) ;
out <- ode(Y, times, func = "derivs", parms = parameters,
           dllname = "testODE", initfunc = "initmod")


plot(out)

Answer 1

编译后的代码不会为 deSolve 在 R 中实现的模型提供不同的结果，除了 atol 限制内的潜在舍入错误和 rtol .

代码中有两个错误的原帖差异的原因。可以按如下方式更正:

1. 声明 static double作为 parms[3];而不是 parms[4]
2. 时间t在 derivs 中是一个指针，即 *t

所以代码读作:

/* file testODE.c */
#include <R.h>
#include <math.h>

static double parms[3];
#define C parms[0] /* left here on purpose */
#define d parms[1]
#define g parms[2]

/* initializer  */
void initmod(void (* odeparms)(int *, double *)) {
  int N=3;
  odeparms(&N, parms);
}

/* Derivatives and 1 output variable */
void derivs (int *neq, double *t, double *y, double *ydot,
             double *yout, int *ip) {
  ydot[0] = -d * exp(-g * *t) * y[0];
}

这里是两个模拟之间的比较，有点适应和概括:

library(deSolve)

testODE <- function(t, y, parameters){
  with(
    as.list(c(y, parameters)),{
      dContamination <- -d * exp(-g * t) * contamination
      return(list(dContamination))
    }
  )
}

system("R CMD SHLIB testODE.c")
dyn.load("testODE.dll")

parameters <- c(c = -8/3, d = -10, g =  28)
Y          <- c(contamination = 1200)
times      <- seq(0, 6, by = 0.01)

out1 <- ode(Y, times, testODE,
            parms = parameters, method = "radau", atol = 1e-4, rtol = 1e-4)
out2 <- ode(Y, times, func = "derivs", dllname = "testODE", initfunc = "initmod",
            parms = parameters, method = "radau", atol = 1e-4, rtol = 1e-4)

plot(out1, out2)          # no visible difference
summary(out1 - out2)      # differences should be (close to) zero

dyn.unload("testODE.dll") # always unload before editing .c file !!

注意:设置 .dll或 .so根据您的操作系统，或使用 .Platform$dynlib.ext 检测它.