1

以下代码代表我的实际问题的简化版本：

# Libraries
library(desolve)
library(parallel)

# End time
max_time <- 300

# Time vector
time <- seq(0, max_time, 0.1)

# Number of events
n_events <- c(6, 60, 600, 6000)

# p1 values
p1_vars <- 10^c(-4, -2, 0, 2)

# Parameters
p <- data.frame(p1 = 1.14*24,
                p2 = 0.14*24,
                p3 = 0.06*24,
                p4 = 0.5,
                p5 = 0.05,
                p6 = 0,
                p7 = 0.218,
                p8 = 0.477,
                p9 = 0.5/6)

# Initial conditions
x0 <- c(x1 = 0,
        x2 = 0,
        x3 = 2e11,
        x4 = 8e11)

# Model
sys <- function(t, y, p) {
  dy = numeric(4)

  dy[1] = p$p2*y[2] - (p$p3 + p$p1)*y[1]
  dy[2] = p$p3*y[1] - p$p2*y[2]
  dy[3] = (p$p4 - p$p7 - p$p8)*y[3] + p$p5*y[4] - p$p9*y[1]*y[3]
  dy[4] = p$p7*y[3] - (p$p5 + p$p6)*y[4]
  list(dy)
}

# Number of corse to parallelise
n_cores <- detectCores() - 1

# Initialise
input.df <- data.frame(eventfreq = rep(n_events/max_time, length(p1_vars)),
                       par = rep(p1_vars, each = length(n_events)))

res.df <- data.frame(eventfreq = NaN,
                     par = NaN,
                     out = NA)

# Parallel function
par_fcn <- function(i) {

  # p1 value
  p$p1 <- input.df$par[i]

  # Event value
  ev_val <- 1000/input.df$eventfreq[i]*max_time

  # Events
  ev <- data.frame(var = 'x1',
                   time = seq(0, max_time, by = 1/input.df$eventfreq[i]),
                   value = 100,
                   method = 'add')

  # Solve ODE
  sim.df <- as.data.frame(ode(y = x0, times = time, func = sys, p = p, events = list(data = ev)))

  # Fill
  res.df$eventfreq <- input.df$eventfreq[i]
  res.df$par <- input.df$par[i]
  res.df$out <- sim.df$x3 + sim.df$x4

  # Return
  return(res.df)
}

# Evaluate in parallel
res_parallel <- mclapply(seq_len(nrow(input.df)), par_fcn, mc.cores = n_cores)
out.df <- do.call(rbind, res_parallel)

如您所见，我尝试使用desolve一系列参数值和不同数量的事件来求解 ODE。尽管我将代码并行化并将其分布在 7 个内核上，但我的实际问题仍然需要很长时间。我已将 ODE 求解器确定为瓶颈，尤其是当事件数量变大时。

有什么想法可以加快此代码的速度吗？

我也一直在研究根据小插图R Package deSolve:Writing Code in Compiled Languages在编译代码中实现这一点。但我对 C 语言太缺乏经验，无法完全靠自己完成。

这是我到目前为止得到的：

/*file model.c */
#include <R.h>
static double p[9];
#define p1 p[0]
#define p2 p[1]
#define p3 p[2]
#define p4 p[3]
#define p5 p[4]
#define p6 p[5]
#define p7 p[6]
#define p8 p[7]
#define p9 p[8]

/* initialiser */
void initmod(void (*odeparms)(int *, double *))
{
    int N = 9;
    odeparms(&N, p);
}

/* Events */
void event(int *n, double *t, double *y)
{
    y[0] = y[0] + 100;
}

/* Derivatives */
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] = p2*y[1] - (p3 + p1)*y[0];
    ydot[1] = p3*y[0] - p2*y[1];
    ydot[2] = (p4 - p7 - p8)*y[2] + p5*y[3] - p9*y[0]*y[2];
    ydot[3] = p7*y[2] - (p5 + p6)*y[3];

    yout[0] = y[0];
    yout[1] = y[2] + y[3];
}

/* END file model.c */

其中，我尤其无法理解如何将n_eventsandp1_vars参数传递给 C 代码。

任何帮助或提示将不胜感激。

4

0 回答 0