Ускорьте параллельное решение событий, содержащих ODE

Question

Следующий код представляет собой упрощенную версию моей актуальной проблемы:

# 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)

Как видите, я пытаюсь решить ODE с desolve для диапазона значений параметров и различного количества событий. Несмотря на то, что я распараллелил код и распределил его по 7 ядрам, моя настоящая проблема все еще занимает очень много времени. В качестве узкого места я определил решатель ODE, особенно когда количество событий становится большим.

Есть ли идеи, как этот код может быть ускорен?

Я также пытался реализовать это в скомпилированном коде в соответствии с виньеткой R Package deSolve: Написание кода на скомпилированных языках. Но я слишком неопытен в Си, чтобы делать это самостоятельно.

Вот что я получил так далеко:

/*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_events а также p1_vars параметр к коду.

Любая помощь или советы будут с благодарностью.