From 357b70cbda678f5d75db8c17a45bfacf7ee0fe5c Mon Sep 17 00:00:00 2001
From: julien arino <julien.arino@gmail.com>
Date: Tue, 28 Mar 2023 12:44:51 +0200
Subject: [PATCH] update

---
 CODE/Julien/SIS_CTMC_parallel.R             | 172 ++++++++++++++++++
 CODE/Julien/SIS_CTMC_parallel_multiple_R0.R | 189 ++++++++++++++++++++
 CODE/Julien/SIS_DTMC_with_markovchain.R     |  55 ++++++
 CODE/Julien/SLIRS_metapop.R                 | 142 +++++++++++++++
 CODE/Julien/simulate_birth_death_CTMC.R     |  85 +++++++++
 5 files changed, 643 insertions(+)
 create mode 100644 CODE/Julien/SIS_CTMC_parallel.R
 create mode 100644 CODE/Julien/SIS_CTMC_parallel_multiple_R0.R
 create mode 100644 CODE/Julien/SIS_DTMC_with_markovchain.R
 create mode 100644 CODE/Julien/SLIRS_metapop.R
 create mode 100644 CODE/Julien/simulate_birth_death_CTMC.R

diff --git a/CODE/Julien/SIS_CTMC_parallel.R b/CODE/Julien/SIS_CTMC_parallel.R
new file mode 100644
index 0000000..62ae46b
--- /dev/null
+++ b/CODE/Julien/SIS_CTMC_parallel.R
@@ -0,0 +1,172 @@
+# Example CTMC simulation of a simple SIS model, parallel version
+library(parallel)
+library(GillespieSSA2)
+library(deSolve)
+
+# Source a file with a few helpful functions for plotting (nice axes labels, crop figure)
+source("useful_functions.R")
+
+# Need a function that runs one simulation and returns a result. While we're at it,
+# we also return an interpolated solution
+run_one_sim = function(params) {
+  IC <- c(S = (params$Pop-params$I_0), I = params$I_0)
+  params_local <- c(gamma = params$gamma, beta = params$beta)
+  reactions <- list(
+    # propensity function effects name for reaction
+    reaction("beta*S*I", c(S=-1,I=+1), "new_infection"),
+    reaction("gamma*I", c(S=+1,I=-1), "recovery")
+  )
+  set.seed(NULL)
+  sol <- ssa(
+    initial_state = IC,
+    reactions = reactions,
+    params = params_local,
+    method = ssa_exact(),
+    final_time = params$t_f,
+    log_firings = TRUE    # This way we keep track of events as well
+  )
+  # Interpolate result (just I will do)
+  wanted_t = seq(from = 0, to = params$t_f, by = 0.01)
+  sol$interp_I = approx(x = sol$time, y = sol$state[,"I"], xout = wanted_t)
+  names(sol$interp_I) = c("time", "I")
+  # Return result
+  return(sol)
+}
+
+# The ODE, in order to compare with solutions to the CTMC
+rhs_SIS_ODE = function(t, x, p) {
+  with(as.list(x), {
+    dS = p$gamma*I-p$beta*S*I
+    dI = p$beta*S*I-p$gamma*I
+    list(c(dS, dI))
+  })
+}
+
+# To run in parallel, it useful to put parameters in a list
+params = list()
+params$Pop = 100
+params$gamma = 1/5
+params$R_0 = 1.5
+params$t_f = 100
+params$I_0 = 1
+# R0 would be (beta/gamma)*S0, so beta=R0*gamma/S0
+params$beta = params$R_0*params$gamma/(params$Pop-params$I_0)
+# Number of simulations. We may want to save all simulation parameters later,
+# so we add it here
+params$number_sims = 50
+
+
+# Detect number of cores (often good to use all but 1, i.e. detectCores()-1)
+# Here, we also test if we are on a 3990X or similar, as R doesn't do 128 threads
+nb_cores <- detectCores()
+if (nb_cores > 124) {
+  nb_cores = 124
+}
+# Initiate cluster
+cl <- makeCluster(nb_cores)
+# Export needed library to cluster
+clusterEvalQ(cl,{
+  library(GillespieSSA2)
+})
+# Export needed variable and function to cluster
+clusterExport(cl,
+              c("params",
+                "run_one_sim"),
+              envir = .GlobalEnv)
+# Run computation
+tictoc::tic()
+SIMS = parLapply(cl = cl, 
+                 X = 1:params$number_sims, 
+                 fun =  function(x) run_one_sim(params))
+tictoc::toc()
+stopCluster(cl)
+
+# The following is if running iteratively rather than in parallel
+if (FALSE) {
+  SIMS = lapply(X = 1:params$number_sims, 
+                FUN =  function(x) run_one_sim(params))
+}
+
+# To compute means of trajectories, we have to work a little. There are two
+# types of means: regular mean, using all trajectories (even those having reached
+# zero) and mean conditioned on non extinction (mean of trajectories not absorbed
+# at zero during the time interval under consideration). We prepare a list that 
+# will contain both means.
+mean_I = list(time = SIMS[[1]]$interp_I$time,
+              I_all = rep(0, length(SIMS[[1]]$interp_I$time)),
+              I_no_extinction = rep(0, length(SIMS[[1]]$interp_I$time)))
+nb_sims_with_NA = 0
+for (i in 1:params$number_sims) {
+  if (any(is.na(SIMS[[i]]$interp_I$I))) {
+    # If we do not condition on extinction, we need to set NAs to 0
+    tmp = SIMS[[i]]$interp_I$I
+    tmp[which(is.na(tmp))] = 0
+    mean_I$I_all = mean_I$I_all+tmp
+    nb_sims_with_NA = nb_sims_with_NA+1
+  } else {
+    mean_I$I_all = mean_I$I_all + SIMS[[i]]$interp_I$I
+    mean_I$I_no_extinction = mean_I$I_no_extinction + SIMS[[i]]$interp_I$I
+  }
+}
+mean_I$I_all = mean_I$I_all/params$number_sims
+mean_I$I_no_extinction = mean_I$I_no_extinction/(params$number_sims-nb_sims_with_NA)
+
+# Now simulate the ODE
+IC <- c(S = (params$Pop-params$I_0), I = params$I_0)
+sol_ODE = ode(y = IC,
+              func = rhs_SIS_ODE,
+              times = seq(from = 0, to = params$t_f, by = 0.1),
+              parms = params)
+
+# Determine maximum value of I for plot
+I_max = max(unlist(lapply(SIMS, function(x) max(x$interp_I$I, na.rm = TRUE))))
+# Prepare y-axis for human readable form
+y_axis = make_y_axis(c(0, I_max))
+
+# Are we plotting for a dark background
+plot_blackBG = TRUE
+if (plot_blackBG) {
+  colour = "white"
+} else {
+  colour = "black"
+}
+
+# Plot
+png(file = "../FIGS/many_CTMC_sims_with_means.png",
+    width = 1200, height = 800, res = 200)
+if (plot_blackBG) {
+  par(bg = 'black', fg = 'white') # set background to black, foreground white
+}
+plot(mean_I$time, SIMS[[1]]$interp_I$I,
+     type = "l", lwd = 0.2, 
+     ylim = c(0, max_I), xaxs = "i",
+     col.axis = colour, cex.axis = 1.25,
+     col.lab = colour, cex.lab = 1.1,
+     yaxt = "n",
+     xlab = "Time (days)", ylab = "Prevalence")
+for (i in 2:params$number_sims) {
+  lines(mean_I$time, SIMS[[i]]$interp_I$I,
+        type = "l", lwd = 0.2)
+}
+lines(mean_I$time, mean_I$I_all,
+      type = "l",
+      lwd = 5, col = "darkorange4")
+lines(mean_I$time, mean_I$I_no_extinction,
+      type = "l",
+      lwd = 5, col = "red")
+lines(sol_ODE[,"time"], sol_ODE[,"I"],
+      type = "l",
+      lwd = 5, col = "dodgerblue4")
+legend("topleft",
+       legend = c("Solutions", "Mean", 
+                  "Mean (not extinct)", "ODE"),
+       cex = 0.6,
+       col = c(ifelse(plot_blackBG, "white", "black"), 
+               "darkorange4", "red", "dodgerblue"),
+       lty = c(1,1,1,1), lwd = c(0.5, 2.5, 2.5, 2.5))
+axis(2, at = y_axis$ticks, labels = y_axis$labels, 
+     las = 1,
+     col.axis = colour,
+     cex.axis = 0.75)
+dev.off()
+crop_figure(file = "../FIGS/many_CTMC_sims_with_means.png")
diff --git a/CODE/Julien/SIS_CTMC_parallel_multiple_R0.R b/CODE/Julien/SIS_CTMC_parallel_multiple_R0.R
new file mode 100644
index 0000000..d5de29e
--- /dev/null
+++ b/CODE/Julien/SIS_CTMC_parallel_multiple_R0.R
@@ -0,0 +1,189 @@
+# Example CTMC simulation of a simple SIS model, parallel version, with multiple R_0 values
+# Here, we just want to know how many simulations see extinction, so we do minimal post-processing
+
+# BEWARE !!!
+#
+# With 100 sims and the number of values used, this can either be very lengthy on a machine with low thread
+# count or very expensive in RAM on a machine with high thread count..
+
+library(GillespieSSA2)
+library(parallel)
+library(dplyr)
+library(latex2exp)
+
+# Source a file with a few helpful functions for plotting (nice axes labels, crop figure)
+source("useful_functions.R")
+
+# Make data frame from list of results
+# Assumes all list entries are similar
+make_df_from_list = function(in_list) {
+  names_list = names(in_list[[1]])
+  tmp = list()
+  for (n in names_list) {
+    tmp[[n]] = unlist(lapply(in_list, function(x) x[[n]]))
+  }
+  OUT = do.call(cbind.data.frame, tmp)
+  return(OUT)
+}
+
+# Need a function that runs one simulation and returns a result. While we're at it,
+# we also return an interpolated solution
+run_one_sim = function(params_vary, params) {
+  R_0 = params_vary$R_0
+  I_0 = params_vary$I_0
+  # Initial conditions
+  IC <- c(S = (params$Pop-I_0), I = I_0)
+  # R0=(beta/gamma)*S0, so beta=R0*gamma/S0
+  beta = R_0*params$gamma/(params$Pop-params$I_0)
+  params_local <- c(gamma = params$gamma, beta = beta)
+  reactions <- list(
+    # propensity function effects name for reaction
+    reaction("beta*S*I", c(S=-1,I=+1), "new_infection"),
+    reaction("gamma*I", c(S=+1,I=-1), "recovery")
+  )
+  set.seed(NULL)
+  sol <- ssa(
+    initial_state = IC,
+    reactions = reactions,
+    params = params_local,
+    method = ssa_exact(),
+    final_time = params$t_f
+  )
+  # If the final time is less than t_f, we've hit extinction
+  if (sol$state[dim(sol$state)[1],"I"] == 0) {
+    sol$extinct = TRUE
+  } else {
+    sol$extinct = FALSE
+  }
+  # To lighten the load, select only a few items to return
+  OUT = list()
+  OUT$R_0 = R_0
+  OUT$I_0 = I_0
+  OUT$final_time = sol$time[length(sol$time)]
+  OUT$final_I = sol$state[length(sol$time),"I"]
+  OUT$extinct = sol$extinct
+  # Return result
+  return(OUT)
+}
+
+# To run in parallel, it useful to put parameters in a list
+params = list()
+params$Pop = 100
+params$gamma = 1/5
+params$R_0 = 1.5
+params$t_f = 100
+params$I_0 = 1
+# R0=(beta/gamma)*S0, so beta=R0*gamma/S0
+params$beta = params$R_0*params$gamma/(params$Pop-params$I_0)
+# Number of simulations. We may want to save all simulation parameters later,
+# so we add it here
+params$number_sims = 500
+
+# Detect number of cores (often good to use all but 1, i.e. detectCores()-1)
+# Here, we also test if we are on a 3990X or similar, as R doesn't do 128 threads
+nb_cores <- detectCores()
+if (nb_cores > 124) {
+  nb_cores = 124
+}
+# Values of I_0 we consider (we do those sequentially)
+values_I_0 = c(1,2,3,5,10)
+# Run main computation loop (iterative part)
+for (I_0 in values_I_0) {
+  writeLines(paste0("I_0=",I_0))
+  # To process efficiently in parallel, we make a list with the different parameter values
+  # we want to change, which will fed by parLapply to run_one_sim
+  params_vary = list()
+  # Vary R_0 and I_0
+  i = 1
+  for (R_0 in c(seq(0.5, 3, by = 0.05))) {
+    for (j in 1:params$number_sims) {
+      params_vary[[i]] = list()
+      params_vary[[i]]$R_0 = R_0
+      params_vary[[i]]$I_0 = I_0
+      i = i+1
+    }
+  }
+  # Initiate cluster
+  cl <- makeCluster(nb_cores)
+  # Export needed library to cluster
+  clusterEvalQ(cl,{
+    library(GillespieSSA2)
+  })
+  # Export needed variable and function to cluster
+  clusterExport(cl,
+                c("params",
+                  "run_one_sim",
+                  "params_vary"),
+                envir = .GlobalEnv)
+  # Run main computation (parallel part)
+  SIMS = parLapply(cl = cl, 
+                 X = params_vary, 
+                 fun =  function(x) run_one_sim(x, params))
+  stopCluster(cl)
+  saveRDS(SIMS, file = sprintf("SIMS_I0_%02d.Rds", I_0))
+}
+
+
+#SIMS = readRDS(file = "SIMS_I0_1.Rds")
+
+# Use dplyr syntax: count the number of extinctions (TRUE and FALSE)
+# after grouping by R_0 value
+results = make_df_from_list(SIMS) %>%
+  count(I_0, R_0, extinct)
+# Take the result and prepare to negate those where n=params$number_sims
+tmp = results %>%
+  filter(n == params$number_sims)
+tmp_insert = tmp
+tmp_insert$extinct = ifelse(tmp$extinct == TRUE, FALSE, TRUE)
+tmp_insert$n = rep(0, dim(tmp_insert)[1])
+# Now inject this result and keep extinct==TRUE
+results = rbind(results, tmp_insert) %>%
+  filter(extinct == TRUE) %>%
+  arrange(I_0, R_0)
+results$pct_extinct = results$n/params$number_sims*100
+
+# Values I_0
+values_I_0 = unique(results$I_0)
+
+# Are we plotting for a dark background
+plot_blackBG = TRUE
+if (plot_blackBG) {
+  colour = "white"
+} else {
+  colour = "black"
+}
+
+# Prepare y-axis for human readable form
+y_axis = make_y_axis(c(0, 100))
+
+# Plot
+png(file = "../FIGS/extinctions_fct_R0_I0.png",
+    width = 1200, height = 800, res = 200)
+if (plot_blackBG) {
+  par(bg = 'black', fg = 'white') # set background to black, foreground white
+}
+tmp = results %>%
+  filter(I_0 == values_I_0[1])
+plot(tmp$R_0, tmp$pct_extinct,
+     type = "b", lwd = 2,
+     ylim = c(0, 100), yaxt = "n",
+     col.axis = colour, cex.axis = 1.25,
+     col.lab = colour, cex.lab = 1.1,
+     xlab = TeX("$R_0$"), ylab = "Percentage extinctions")
+for (i in 2:length(values_I_0)) {
+  tmp = results %>%
+    filter(I_0 == values_I_0[i])
+  lines(tmp$R_0, tmp$pct_extinct,
+        type = "l", lwd = 2, lty = i)
+}
+abline(v = 1, lty = "dotted")
+# legend("topright",
+#        legend = TeX(sprintf("$I_0$=%d", values_I_0)),
+#        lty = 1:length(values_I_0), lwd = rep(2, length(values_I_0)))
+axis(2, at = y_axis$ticks, labels = y_axis$labels, 
+     las = 1,
+     col.axis = colour,
+     cex.axis = 0.75)
+dev.off()
+crop_figure(file = "../FIGS/extinctions_fct_R0_I0.png")
+
diff --git a/CODE/Julien/SIS_DTMC_with_markovchain.R b/CODE/Julien/SIS_DTMC_with_markovchain.R
new file mode 100644
index 0000000..6049686
--- /dev/null
+++ b/CODE/Julien/SIS_DTMC_with_markovchain.R
@@ -0,0 +1,55 @@
+# Example simulation of a simple SIS model
+library(markovchain)
+
+# Source a file with a few helpful functions for plotting (nice axes labels, crop figure)
+source("useful_functions.R")
+
+# Total population
+Pop = 5
+# Initial number of infectious
+I_0 = 2
+# Parameters
+gamma = 1/5
+R_0 = 1.5
+# R0 would be (beta/gamma)*S0, so beta=R0*gamma/S0
+beta = R_0*gamma/(Pop-I_0)
+# Time step
+Delta_t = 1
+
+# Make the transition matrix
+T = mat.or.vec(nr = (Pop+1), nc = (Pop+1))
+for (row in 2:Pop) {
+  I = row-1
+  mv_right = gamma*I*Delta_t # Recoveries
+  mv_left = beta*I*(Pop-I)*Delta_t # Infections
+  T[row,(row-1)] = mv_right
+  T[row,(row+1)] = mv_left
+}
+# Last row only has move left
+T[(Pop+1),Pop] = gamma*(Pop)*Delta_t
+# Check that we don't have too large values
+if (max(rowSums(T))>1) {
+  T = T/max(rowSums(T))
+}
+diag(T) = 1-rowSums(T)
+
+
+mcSIS <- new("markovchain", 
+             states = sprintf("I_%d", 0:Pop),
+             transitionMatrix = T,
+             name = "SIS")
+
+# Show some info about the chain
+summary(mcSIS)
+
+# The vector of steady states. Here, all mass should be in I_0 (the first state)
+steadyStates(mcSIS)
+
+# Here, we only have one absorbing state, so this is trivial
+meanRecurrenceTime(mcSIS)
+
+# Hitting probabilities
+hittingProbabilities(mcSIS)
+
+# Mean absorption time
+meanAbsorptionTime(mcSIS)
diff --git a/CODE/Julien/SLIRS_metapop.R b/CODE/Julien/SLIRS_metapop.R
new file mode 100644
index 0000000..128adb6
--- /dev/null
+++ b/CODE/Julien/SLIRS_metapop.R
@@ -0,0 +1,142 @@
+# SLIRS_metapop
+#
+# R script that runs integration of the sample p-SLIRS metapopulation
+# model in the paper by Arino (2017).
+#
+# Julien Arino
+# February 2017
+
+# Note that without further adaptation, the results provided by this
+# script can be expected to be completely unrealistic, the intent is only 
+# to provide a workable canvas.
+
+
+# Load required library. 
+# You may need to install this library using install.packages("deSolve").
+library(deSolve)
+
+
+# The vector field
+SLIRS_metapop_rhs <- function (t, x, p) {
+  # We don't use with(as.list(x)) because we need to name the vectors anyway
+  S = x[p$idx_S]
+  L = x[p$idx_L]
+  I = x[p$idx_I]
+  R = x[p$idx_R]
+  # We need N as well for the incidence function
+  N = S+L+I+R
+  # For simplicity, we pre-fill the vector of incidence functions.
+  Phi = p$beta*S*I/(N^p$eta)
+  # Now we set the values of the derivatives. 
+  dS = p$b+p$nu*R-p$d*S-Phi+p$MS%*%S
+  dL = Phi-(p$epsilon+p$d)*L+p$ML%*%L
+  dI = p$epsilon*L-(p$gamma+p$d)*I+p$MI%*%I
+  dR = p$gamma*I-(p$nu+p$d)*R+p$MR%*%R
+  dx = list(c(dS, dL, dI, dR))
+  return(dx)
+}
+
+# Function for computing R0 given we know the 21 entry in V inverse
+R0_fct = function(p) {
+  Sstar = solve(diag(p$d)-p$MS) %*% p$b
+  # Case of mass action incidence
+  #F12 = diag(as.vector(p$beta*Sstar))
+  # Case of proportional incidence
+  F12 = diag(as.vector(p$beta*(Sstar^(1-p$eta))))
+  Vinv21 = solve(diag(p$epsilon+p$d)-p$ML) %*% 
+    diag(p$epsilon) %*% solve(diag(p$gamma+p$d)-p$MI)
+  FVinv = F12 %*% Vinv21
+  ev = eigen(FVinv)
+  Rzero = max(abs(ev$values))
+  return(Rzero)
+}
+
+# The following data originates from Bluedot.global and shows the average
+# daily number of air passengers flying between each of the countries
+# listed below. Average computed over one (unspecified) month.
+# pop contains the population of the countries.
+pop = c(34.017,1348.932,1224.614,173.593,93.261)*1e6
+countries = c("Canada","China","India","Pakistan","Philippines")
+T = matrix(data = c(0,1268,900,489,200,
+                    1274,0,678,859,150,
+                    985,703,0,148,58,
+                    515,893,144,0,9,
+                    209,174,90,2,0),
+           nrow =5, ncol = 5, byrow = TRUE)
+
+# We use a list to store parameters.
+p = list()
+
+# In order to compute the movement rates, we use the approximation
+# explained in Arino & Portet, JMB 2015.
+p$M = mat.or.vec(nr = dim(T)[1], nc = dim(T)[2])
+for (from in 1:5) {
+  for (to in 1:5) {
+    p$M[to,from] = -log(1-T[from,to]/pop[from])
+  }
+  p$M[from,from] = 0
+}
+p$MS = p$M-diag(colSums(p$M))
+# Assume movement is equal for all compartments.
+p$ML = p$MS
+p$MI = p$MS
+p$MR = p$MS
+# Number of patches. (MS should be square.)
+p$P = dim(p$MS)[1]
+
+# Set other parameters. Values are chosen mostly arbitrarily.
+p$beta = rep(0.1,p$P)
+p$nu = rep(0.1,p$P)
+p$epsilon = rep(0.01,p$P)
+p$d = rep(1/(70*365.25),p$P) # Could be taken based on World Bank data to be more accurate
+p$gamma = rep(0.05,p$P)
+p$eta = rep(1,p$P)
+# We use the numerical trick explained in Arino & Portet, JMB 2015, in
+# order to work with a constant total population. Note that this can lead
+# to negative birth rates and should only be used if one intends to start
+# with an initial condition with the total population equal to its
+# equilibrium value.
+p$b = (diag(p$d)-p$MS) %*% pop
+
+# Save index of state variable types in state variables vector (we have to use 
+# a vector and thus, for instance, the name "S" needs to be defined)
+p$idx_S = 1:p$P
+p$idx_L = (p$P+1):(2*p$P)
+p$idx_I = (2*p$P+1):(3*p$P)
+p$idx_R = (3*p$P+1):(4*p$P)
+# Set initial conditions. For example, we start with 2 infectious
+# individuals in Canada.
+L0 = mat.or.vec(p$P,1)
+I0 = mat.or.vec(p$P,1)
+I0[1] = 2
+R0 = mat.or.vec(p$P,1)
+S0 = pop-(L0+I0+R0)
+# Vector of initial conditions to be passed to ODE solver.
+IC = c(S = S0, L = L0, I = I0, R = R0)
+# Time span of the simulation (5 years here).
+tspan = seq(from = 0, to = 5*365.25, by = 0.1)
+
+# Call of the ODE solver.
+sol <- ode(y = IC, times = tspan, func = SLIRS_metapop_rhs, parms = p)
+
+# Put solution in an easier to use form.
+times = sol[,"time"]
+S = sol[,p$idx_S]
+L = sol[,p$idx_L]
+I = sol[,p$idx_I]
+R = sol[,p$idx_R]
+N = S+L+I+R
+
+# A sample (and simple) plot: number infected per 100,000 inhabitants.
+xlim = range(times)
+ylim = c(0,max(I/N*1e5))
+plot(0,
+     xlim = xlim, ylim = ylim,
+     xlab = "Time (days)", ylab = "Number infectious per 100,000")
+for (i in 1:p$P) {
+  lines(times,I[,i]/N[,i]*1e5, col = i)
+}
+legend("topleft",legend = countries,col = 1:p$P, lty = 1)
+
+# Finally, compute the value of R0
+R0_fct(p)
diff --git a/CODE/Julien/simulate_birth_death_CTMC.R b/CODE/Julien/simulate_birth_death_CTMC.R
new file mode 100644
index 0000000..b14150a
--- /dev/null
+++ b/CODE/Julien/simulate_birth_death_CTMC.R
@@ -0,0 +1,85 @@
+# Running a birth-death process as a CTMC using classic Gillespie
+# 
+# This demos both the classic Gillespie algorithm and how things can go bad
+# because of shrinking inter-event times.
+
+b = 0.025   # Birth rate
+d = 0.01   # Death rate
+t_0 = 0    # Initial time
+N_0 = 100  # Initial population
+
+# Vectors to store time and state. Initialise with initial condition.
+t = t_0
+N = N_0
+
+t_f = 1000  # Final time
+
+# We'll track the current time and state (could also just check last entry in t
+# and N, but will take more operations)
+t_curr = t_0
+N_curr = N_0
+
+while (t_curr<=t_f) {
+  xi_t = (b+d)*N_curr
+  # The exponential number generator does not like a rate of 0 (when the 
+  # population crashes), so we check if we need to quit
+  if (N_curr == 0) {
+    break
+  }
+  tau_t = rexp(1, rate = xi_t)
+  t_curr = t_curr+tau_t
+  v = c(b*N_curr, xi_t)/xi_t
+  zeta_t = runif(n = 1)
+  pos = findInterval(zeta_t, v)+1
+  switch(pos,
+         { 
+           # Birth
+           N_curr = N_curr+1
+         },
+         {
+           # Death
+           N_curr = N_curr-1
+         })
+  N = c(N, N_curr)
+  t = c(t, t_curr)
+}
+
+# Make a title for all figures
+main_title = paste0("b=",b,", d=",d)
+fig_title = gsub("\\.", "_", paste0("b=",b,"__d=",d))
+
+# Do we output the results?
+OUTPUT_PLOT = TRUE
+
+# Plot the result
+if (OUTPUT_PLOT) {
+  png(filename = paste0("../FIGS/CTMC_birth_death_sol_", fig_title, ".png"),
+      width = 800, height = 600, res = 120)
+}
+plot(t, N, type = "l",
+     main = main_title)
+abline(h = N_0, lty = 2)
+if (OUTPUT_PLOT) {
+  dev.off()
+  crop_figure(paste0("../FIGS/CTMC_birth_death_sol_", fig_title, ".png"))
+}
+
+# Plot the inter-event time. Diff removes one entry, beware..
+plot(1:(length(t)-1), diff(t),
+     xlab = "Jump #", ylab = "Inter-event time",
+     main = main_title,
+     type = "l")
+
+# Plot the inter-event time. Diff removes one entry, beware..
+if (OUTPUT_PLOT) {
+  png(filename = paste0("../FIGS/CTMC_birth_death_ie_vs_t_", fig_title, ".png"),
+      width = 800, height = 600, res = 120)
+}
+plot(t[2:length(t)], diff(t),
+     xlab = "Time", ylab = "Inter-event time",
+     main = main_title,
+     type = "l")
+if (OUTPUT_PLOT) {
+  dev.off()
+  crop_figure(paste0("../FIGS/CTMC_birth_death_ie_vs_t_", fig_title, ".png"))
+}