title: Speed up parallelised solving of ODE containing events
body: <p>The following code represents a simplified version of my actual problem:</p>
  
  <pre class="lang-r prettyprint-override"><code># Libraries
  library(desolve)
  library(parallel)
  
  # End time
  max_time &lt;- 300
  
  # Time vector
  time &lt;- seq(0, max_time, 0.1)
  
  # Number of events
  n_events &lt;- c(6, 60, 600, 6000)
  
  # p1 values
  p1_vars &lt;- 10^c(-4, -2, 0, 2)
  
  # Parameters
  p &lt;- 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 &lt;- c(x1 = 0,
          x2 = 0,
          x3 = 2e11,
          x4 = 8e11)
  
  # Model
  sys &lt;- 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 &lt;- detectCores() - 1
  
  # Initialise
  input.df &lt;- data.frame(eventfreq = rep(n_events/max_time, length(p1_vars)),
                         par = rep(p1_vars, each = length(n_events)))
  
  res.df &lt;- data.frame(eventfreq = NaN,
                       par = NaN,
                       out = NA)
  
  # Parallel function
  par_fcn &lt;- function(i) {
  
    # p1 value
    p$p1 &lt;- input.df$par[i]
  
    # Event value
    ev_val &lt;- 1000/input.df$eventfreq[i]*max_time
  
    # Events
    ev &lt;- data.frame(var = 'x1',
                     time = seq(0, max_time, by = 1/input.df$eventfreq[i]),
                     value = 100,
                     method = 'add')
  
    # Solve ODE
    sim.df &lt;- as.data.frame(ode(y = x0, times = time, func = sys, p = p, events = list(data = ev)))
  
    # Fill
    res.df$eventfreq &lt;- input.df$eventfreq[i]
    res.df$par &lt;- input.df$par[i]
    res.df$out &lt;- sim.df$x3 + sim.df$x4
  
    # Return
    return(res.df)
  }
  
  # Evaluate in parallel
  res_parallel &lt;- mclapply(seq_len(nrow(input.df)), par_fcn, mc.cores = n_cores)
  out.df &lt;- do.call(rbind, res_parallel)
  </code></pre>
  
  <p>As you can see, I try to solve an ODE with <code>desolve</code> for a range of parameter values and different number of events. Even though, I parallelised the code and spread it across 7 cores, my actual problem still takes very long. I've identified the ODE solver as the bottleneck, especially when the number of events gets large.</p>
  
  <p>Are there any ideas how this code could be speeded up?</p>
  
  <p>I've also been looking into implementing this in compiled code according to the vignette <a href="https://cran.r-project.org/web/packages/deSolve/vignettes/compiledCode.pdf" rel="nofollow noreferrer">R Package deSolve: Writing Code in Compiled Languages</a>. But I'm too inexperienced in C to do it completely on my own.</p>
  
  <p>Here's what I got so far:</p>
  
  <pre class="lang-c prettyprint-override"><code>/*file model.c */
  #include &lt;R.h&gt;
  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(&amp;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] &lt; 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 */
  </code></pre>
  
  <p>Therein, I especially can't follow the vignette of how to pass the <code>n_events</code> and <code>p1_vars</code> parameter to the C code.</p>
  
  <p>Any help or hints would be greatly appreciated.</p>