README.Rmd

---
title: "ssar: Stochastic Simulation Algorithm in R"
output: github_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

The R package _ssar_ is a fast implementation of Gillespie's Stochastic Simulation Algorithm. It combines R's graphical and statistical capabilities with the speed of C++. In addition, the package allows for simulation of stochastic processes with time-dependent propensity functions Thus, _ssar_ represents an improvement over the previous package available at CRAN [_GillespieSSA_](https://cran.r-project.org/web/packages/GillespieSSA/index.html).

## Is this package for me?

This package is for you if: 

  * You want to simulate Continuous Time Markov Chains (CTMC), Stochastic Compartmental Models (like the ones in chemistry, ecology, epidemiology).
  
  * You want to use Gillespie's Stochastic Simulation Algorithm with time-dependent parameters, random parameters and/or time-dependent random parameters. (Areas might include Bayesian MonteCarlo parameter estimation of some Stochastic Processes or Inhomogeneous Continuous Time Markov Chains)
  
  * You are tired of the current packages being too slow.
  
  * You want to have fun simulating stuff!

## Installation

_ssar_  is still at its developmental stage. You need to install from github:

```{r install, eval  = FALSE}
install.packages("devtools")
devtools::install_github("INSP-RH/ssar")
```

To adecquately run ``ssar``you need to install a `C++`compiler: 

+ On Windows, install [Rtools](https://cran.r-project.org/bin/windows/Rtools/).
+ On Mac, install Xcode from the App store.
+ On Linux, ``sudo apt-get install r-base-dev`` or similar.

## Examples 

After installing you need to call the _ssar_ package for using it. 

```{r library}
library(ssar)
```

### 1. Logistic Growth

First, we set the seed for the simulation. 
```{r, setseed}
set.seed(123)
```


Initial data must be inputed as a matrix. 
```{r Xlogistic}
X          <- matrix(c(N=500), nrow = 1)
```

The propensity vector should also be in matrix form:
```{r vlogistic}
v          <- matrix( c(+1, -1), ncol = 2)
```

The propensity scores must also be a matrix-valued function depdendent on 3 parameters: time (`t`), the state of the system (`X`) and additional parameters (`params`) which we discuss later.
```{r pfunlogistic}
pfun       <- function(t,X,params){ cbind(2 * X[, 1], (1 + 1*X[, 1]/1000)*X[, 1]) }
```


The model runs automatically from _0_ to _1_ conducting _10_ simulations and generating a plot.
```{r simulation1, results = 'hide'}
simulation <- ssa(X, pfun, v)
```

The `nsim` variable specifies the number of simulations. Suppose we want _20_ of them:
```{r simulation2, results = 'hide'}
simulation <- ssa(X, pfun, v, nsim = 20)
```

The `tmin` and `tmax` variables specify the initial time and final time of the process. Suppose we want to simulate from time _2_ to time _10_ with _20_ simulations:

```{r simulation3, results = 'hide'}
simulation <- ssa(X, pfun, v, tmin = 2, tmax = 10, nsim = 20)
```

Plot characteristics can be specified by `title`, `xlab` and `ylab`:
```{r simulation4, results = 'hide'}
simulation <- ssa(X, pfun, v, tmin = 2, tmax = 10, nsim = 20, 
                  title = "Logistic Growth: Example 1", 
                  xlab = "Time", ylab = "Individuals")
```

Making plots can really slow down the process. The option: `plot.sim` when set to `FALSE` allows us to keep the data without making any plot: 
```{r simulation5, eval = FALSE, results = 'hide'}
simulation <- ssa(X, pfun, v, tmin = 2, tmax = 10, nsim = 20, plot.sim = FALSE)
```

The `simulation` dataframe looks like this: 

```{r simulation6, echo = FALSE}
simulation <- ssa(X, pfun, v, tmin = 2, tmax = 10, nsim = 20, plot.sim = FALSE)
head(simulation)
```

### 2. Time-dependent Logistic Growth

Suppose we are using almost the same model as in the previous example:

```{r, setseed2}
set.seed(322)
X          <- matrix(c(N=500), nrow = 1)   #Initial values
v          <- matrix( c(+1, -1), ncol = 2) #Propensity scores
```

But the propensity function now depends on time:
```{r, pfun2}
pfun       <- function(t,X,params){ cbind(2 * X[, 1], 
                                          (2 + sin(t*pi)*X[, 1]/1000)*X[, 1]) }
```

Simulation is done in exactly the same manner as previously done. No change needed!
```{r simulation7, results = 'hide'}
simulation <- ssa(X, pfun, v, tmin = 2, tmax = 10, nsim = 20, 
                  title = "Time-dependent Logistic Growth: Example2", 
                  xlab = "Time", ylab = "Individuals")
```

### 3. Exponential model

This is a new model given by the following parameters: 
```{r, pfun3}
 #Start the parameters
 X          <- matrix(c(N = 10), nrow = 1)
 pfun       <- function(t, X, params){ cbind(1.1 + sin(pi*t/0.01))*X[,1] }
 v          <- matrix( c(+1), ncol=1)
```

```{r, simulation8, results = 'hide'}
 simulation <- ssa(X, pfun, v,
                   title = "Example 3", xlab = "Time", ylab = "Value")
```

The option `maxiter` establishes the maximum number of iterations done by the model before stopping. For example if we wish to know where the model is after _100_ changes set `maxiter = 100`:

```{r, simulation9, results = 'hide'}
 simulation <- ssa(X, pfun, v, maxiter = 100, 
                   title ="Model after 100 changes: Example 3",
                   xlab = "Time", ylab = "Value")
```

The option `print.time` prints to screen at what in time of the simulation we are. For example if the model goes from `tmin = 0` to `tmax = 1` setting `print.time = TRUE` will print at which moment in time the model is simulating: 

```{r, simulation10, results = 'hide', eval = FALSE}
 simulation <- ssa(X, pfun, v, maxiter = 100, tmin = 0, tmax = 1,
                   plot.sim = FALSE, print.time = TRUE)
```

```{r, simulation11, echo = FALSE}
 simulation <- ssa(X, pfun, v, maxiter = 5, tmin = 0, tmax = 1, 
                   plot.sim = FALSE, print.time = TRUE)
```

The option `maxtime` establishes how much computer-time (in seconds) will be used for the model. This is specially useful for models which might take a lot of time to run. In the following example, we run the model for _2_ seconds:

```{r, simulation12, results = 'hide'}
 simulation <- ssa(X, pfun, v, maxtime = 2, 
                   title ="Model after 2 seconds: Example 3",
                   xlab = "Time", ylab = "Value" )
```

### 4. Lotka-Volterra
<a name="lvolterra"></a>
We find it easier to assign the _parameters_ (constants) used by the propensity function as a separate vector. This is done in the following simulation: 

```{r simulation13, results = 'hide'}
#Set seed
set.seed(3289650)

#Get initial parameters
params     <- c(a = 3, b = 0.01, c = 2)
X          <- matrix(c(100, 100), ncol = 2)

#Propensity function
pfun       <- function(t, X, params){ cbind(params[1]*t*X[,1] + 1, 
                                              params[2]*X[,1]*X[,2], 
                                                params[3]*X[,2]) }
#Propensity score
v          <- matrix(c(+1,-1,0,0,+1,-1),nrow=2,byrow=TRUE)

#Simulate
simulation <- ssa(X, pfun, v, params, 
                  title = "Example 4: Time-dependent Lotka-Volterra",
                  xlab = "Time", ylab = "Number of individuals")
```

<a name="keepfilesec"></a> The `ssa` function works by creating a file called *"Temporary\_File\_ssa.txt"*. Setting to `TRUE` option `keep.file` does not remove the temporary file. Furthermore, the option `fname` allows you to rename the file. This option is really helpful if you want to keep a database of your simulation: 

```{r simulation14, results = 'hide'}
simulation <- ssa(X, pfun, v, params, keep.file = TRUE, fname ="My_simulation.txt",
                  plot.sim = FALSE)
```

You can read the file with the `read.table` function: 
```{r simulation15, results = 'hide', eval = FALSE}
sim        <- read.table("My_simulation.txt",  header = TRUE)
```

```{r simulation16, echo = FALSE}
sim        <- read.table("My_simulation.txt", header = TRUE)
head(sim)
```

If you are a `ggplot2` kind of person you can plot easily your simulations:
```{r simulation17}
library(ggplot2)
ggplot(data = sim, aes(x = Time, group = as.factor(Simulation))) +
    geom_line(aes(y = Var1, color = "Prey")) +
    geom_line(aes(y = Var2, color = "Predator")) +
    ggtitle("Example 4: Lotka Volterra with ggplot2") + 
    xlab("Time") + ylab("Individuals") +
    scale_color_manual("Creature", 
                       values = c("Prey" = "deepskyblue4","Predator" = "tomato3"))
```

### 5. Lotka-Volterra with random time-dependent parameters

This is almost the same Lotka-Volterra model; however in this case the parameters _a_ and _b_ are random variables. 


```{r simulation18, results = 'hide'}
#Set seed
set.seed(3289650)

#Get initial parameters
params     <- c(amu = 3, asd = 0.01, bmin = 0.001, bmax = 0.015, c = 2)
X          <- matrix(c(100, 100), ncol = 2)

#Propensity function
pfun       <- function(t, X, params){ cbind(rnorm(1,params[1], params[2])*X[,1] + 1, 
                                            runif(1,params[3],params[4])*X[,1]*X[,2], 
                                                params[5]*X[,2]) }
#Propensity score
v          <- matrix(c(+1,-1,0,0,+1,-1),nrow=2,byrow=TRUE)

#Simulate
simulation <- ssa(X, pfun, v, params, 
                  title = "Example 5: Lotka-Volterra with random variables",
                  xlab = "Time", ylab = "Number of individuals")
```

Notice that the random variables in the model can also be time-dependent: 
```{r simulation19, results = 'hide'}
#Propensity function
pfun       <- function(t, X, params){ 
              cbind(rnorm(1,t + params[1], params[2])*X[,1] + 1, 
                    runif(1,params[3],params[4])*X[,1]*X[,2], params[5]*X[,2]) }

#Simulate
simulation <- ssa(X, pfun, v, params, 
              title = "Example 5: Lotka-Volterra with time-dependent random variables",
              xlab = "Time", ylab = "Number of individuals")
```

### 6. Additional tips for running faster and/or with less memory

Sometimes your model might take a lot of time to run. The following list of options might help you speed it up: 

+ Do not print the current time: `print.time = FALSE`

The fastest way to speed up your code is via the  `file.only` and `kthsave` options

### The `file.only` option

As we said in the [previous section](#keepfilesec), the program generates a Temporary File. The `file.only` option generates the file but does not return any values to `R`nor does it generate a plot. It is meant for making fast simulations in which the user might not be interested in generating a plot inside the function. 

As an example, consider the [Lotka-Volterra model](#lvolterra). 

```{r LVM, results = 'hide'}
#Set seed
set.seed(3289650)

#Get initial parameters
params     <- c(a = 3, b = 0.01, c = 2)
X          <- matrix(c(100, 100), ncol = 2)

#Propensity function
pfun       <- function(t, X, params){ cbind(params[1]*t*X[,1] + 1, 
                                              params[2]*X[,1]*X[,2], 
                                                params[3]*X[,2]) }
#Propensity score
v          <- matrix(c(+1,-1,0,0,+1,-1),nrow=2,byrow=TRUE)

```


Without the `file.only` option:
```{r LVM1, echo = FALSE}
init1   <- Sys.time()
```

```{r LVM2, results = 'hide'}
#Simulate
simulation <- ssa(X, pfun, v, params, plot.sim = FALSE)
```

```{r LVM3, echo = FALSE}
finit1  <- Sys.time()
```

With the `file.only` option:
```{r LVM4, echo = FALSE}
init2   <- Sys.time()
```

```{r LVM5, results = 'hide'}
#Simulate
simulation <- ssa(X, pfun, v, params, file.only = TRUE)
```

```{r LVM6, echo = FALSE}
finit2  <- Sys.time()
```

```{r LVM7, echo = FALSE}
print("********** OVERALL TIME EVALUATION **********")
print(paste0("Normal eval:     ", finit1 - init1))
print(paste0("file.only = TRUE: ", finit2 - init2))
print("*********************************************")
```

This might not look as fast; However in bigger files, it is really important. Additional benchmarks are provided in the [Benchmarking section](#benchmarking)

### The `kthsave` option

The Stochastic Simulation Algorithm computes and saves every transition made in the model. This might not be a problem for short simulations; but in the long run generates large databases which are pretty memory intensive. The `kthsave` option is here to help. 

Consider the following model which is a variant of the SIS model for epidemics: 

```{r simulation20, results = 'hide'}
 #Initial parameters
 k          <-  24576.5529836797
 delta      <-  0.0591113454895868 + 0.208953907151055
 gamma_ct   <-  0.391237630231631
 params     <- c(k = k, delta = delta, gamma_ct = gamma_ct)
 X          <- matrix(c(S = 1000000000, I = 1000), ncol = 2)
 pfun       <- function(t, X, params){

   #Value to return
   matreturn  <- matrix(NA, nrow = length(t), ncol = 6)

   #Create birth function
   lambda     <- function(t){ return(4.328e-4 - (2.538e-7)*t - 
                              (3.189e-7)*sin(2 * t * pi/52) - 
                              (3.812e-7)*cos(2 * t * pi/52) ) }

   #Create death function
   mu         <- function(t){ return(9.683e-5 + (1.828e-8)*t + 
                              (2.095e-6)*sin(2 * t * pi/52) - 
                              (8.749e-6)*cos(2 * t * pi/52))}

   #Create infectives function
   beta_fun   <- function(t){ return( 0.479120824267286 + 
                              0.423263042762498*sin(-2.82494252560096 + 2*t*pi/52) )}

   #Estimate values
   matreturn[,1] <- lambda(t)*(X[,1] + X[,2])
   matreturn[,2] <- mu(t)*X[,1]
   matreturn[,3] <- beta_fun(t)*X[,1]*X[,2]/(1 + params[1]*X[,2])
   matreturn[,4] <- mu(t)*X[,2]
   matreturn[,5] <- params[2]*X[,2]
   matreturn[,6] <- params[3]*X[,2]

   #Return
   return(matreturn)

 }
 v          <- matrix(c(1,-1, -1, 0, 0, 1, 0, 0, 1, -1, -1, -1), nrow = 2, byrow = TRUE)
 tmin       <- 0
 tmax       <- 2
 nsim       <- 100

```

Running _100_ simulations for _2_ days generates over 4 GB of information:

**DO NOT RUN: MIGHT TAKE SEVERAL MINUTES** 
```{r simulation21, eval = FALSE}
 #DO NOT RUN 
 simulation <- ssa(X, pfun, v, params, tmin, tmax, nsim = nsim, print.time = FALSE, 
                    plot.sim = FALSE, keep.file = TRUE)
 #DO NOT RUN 
```

Running the simulation for _52_ days generates over 30 GB of information. In order to speed the program and reduce the simulation time we can only save every _kth_ iteration. The command _kthsave_ does te trick. 

The first 1000 iterations of the model look like this:
```{r time1, echo = FALSE}
init1  <- Sys.time()
```

```{r simulation22, results = 'hide'}
 set.seed(123)
 simulation1 <- ssa(X, pfun, v, params, tmin, tmax, nsim = 10, print.time = FALSE, 
                    plot.sim = FALSE, maxiter = 5000, keep.file = TRUE,
                    fname = "sim1.txt")
 
```

```{r time2, echo = FALSE}
finit1 <- Sys.time()
```

We now consider saving only every _10_ iterations of the model:

```{r time3, echo = FALSE}
init2 <- Sys.time()
```

```{r simulation23, results = 'hide'}
 set.seed(123)
 simulation2 <- ssa(X, pfun, v, params, tmin, tmax, nsim = 10, print.time = FALSE, 
                    plot.sim = FALSE, maxiter = 5000, kthsave = 10, keep.file = TRUE,
                    fname = "sim2.txt")
```

```{r time4, echo = FALSE}
finit2 <- Sys.time()
```

There are almost no noticable differences between the models: 

```{r plotsim23, results = 'hide'}
 ggplot(simulation1, aes(x = Time, y = Var2, group=as.factor(Simulation))) +
    geom_point(data = simulation2, 
               aes(color = "Every 10 values")) + 
    geom_step(data = simulation1, 
               aes(color = "All values"), size = 0.5) + 
    theme(legend.position="none") + theme_bw() +
    ggtitle(paste0("SIS example; Infected cases ", 10, " simulations")) + 
    xlab("Time") + ylab("Individuals") 
```

Changing `kthsave` to _10_ reduces _10_ times the file size. In addition, it almost halves the modeling speed:

```{r time5, echo = FALSE}
print("********** OVERALL TIME EVALUATION **********")
print(paste0("All values:      ", finit1 - init1))
print(paste0("Every 10 values: ", finit2 - init2))
print("*********************************************")
```


## Benchmarking
<a name="benchmarking"></a>
In order to show the advantage or this package over the existing [_GillespieSSA_](https://cran.r-project.org/web/packages/GillespieSSA/index.html) we show several benchmarks. Running this in your computer requires installation of the _GillespieSSA_ package and the _microbenchmark_ package. 

```{r benchmark0, results = 'hide'}
library(microbenchmark)
```

First, we run the program from _GillespieSSA_:

```{r benchmark1}
#Running the program from GillespieSSA
set.seed(1)
parms    <- c(c=0.5)
  x0     <- c(X=10000)
  a      <- c("c*X")
  nu     <- matrix(-1)
  
gilltime <- microbenchmark(
                out1 <- GillespieSSA::ssa(x0,a,nu,parms,tf = 100)
)

```

Notice that running this program results in only one simulation: 
```{r benchmark2, echo = FALSE}
GillespieSSA::ssa.plot(out1)
```

In the case of our model: we can make 5 simulations faster than 1 simulation from _GillespieSSA_:

```{r benchmark3, results = 'hide'}
set.seed(1)

parms    <- c(0.5)
  x0     <- matrix(c(X=10000), ncol = 1)
  pfun   <- function(t,X,params){ return(as.matrix(params[1]*X[,1])) }
  nu     <- matrix(-1)

#Keeping all the information
alltime1  <- microbenchmark(
              out2 <- ssar::ssa(x0, pfun, nu, parms, tmin = 0, 
                                tmax = 100, nsim = 5, plot.sim = FALSE))

#All the iterations and file.only option
alltime2  <- microbenchmark(
              out2 <- ssar::ssa(x0, pfun, nu, parms, tmin = 0, 
                                tmax = 100, nsim = 5, file.only = TRUE))

#Keeping every 10 iterations. This is really fast in comparison.
tentime1  <- microbenchmark(
              out2 <- ssar::ssa(x0, pfun, nu, parms, tmin = 0, 
                                tmax = 100, nsim = 5, plot.sim = FALSE,
                                kthsave = 10))

#10 iterations and file.only option
tentime2  <- microbenchmark(
              out2 <- ssar::ssa(x0, pfun, nu, parms, tmin = 0, 
                                tmax = 100, nsim = 5, file.only = TRUE,
                                kthsave = 10))

```

Notice that _ssar_ creates 5 simulations: 
```{r benchmark5, echo = FALSE, results = 'hide'}
val<-ssar::ssa(x0, pfun, nu, parms, tmin = 0, tmax = 100, nsim = 5, plot.sim = T)
```

The overall times (seconds):
```{r overalltimes, echo = FALSE}
timedata <- as.data.frame(rbind(summary(gilltime[[2]]), 
                                summary(alltime1[[2]]), 
                                summary(alltime2[[2]]),
                                summary(tentime1[[2]]), 
                                summary(tentime2[[2]])))/1.e9
rownames(timedata) <- c("GillespieSSA","All times","All times file.only",
                        "kthsave 10","kthsave 10 and file.only")
colnames(timedata) <- names(summary(alltime1[[2]]))

knitr::kable(timedata, format = "markdown")

```

## Common errors and their meaning

### `pfun` needs to be a matrix valued function

The function `pfun`is not returning a matrix. You can use `as.matrix` or `cbind` (depending on your function) to return a matrix value. As an example: 


```{r error1, eval = FALSE}
#THIS IS INCORRECT (USING c NOT cbind):
pfun       <- function(t,X,params){ c(params[1] *(1 + sin(t))* X[,1], 
                          (params[2] + (params[1]-params[2])*X[,1]/params[3])*X[,1]) }

#THIS IS CORRECT:
pfun       <- function(t,X,params){ cbind(params[1] *(1 + sin(t))* X[,1], 
                          (params[2] + (params[1]-params[2])*X[,1]/params[3])*X[,1] ) }
```

### `xinit` needs to be a matrix

The value `xinit` is probably a vector and not a matrix. 

```{r error2, eval = FALSE}
#THIS IS INCORRECT (USING c NOT cbind):
xinit <- c(X = 1, Y = 2)

#THIS IS CORRECT:
xinit <- matrix(c(X = 1, Y = 2), ncol = 1)
```

### Error in `pfun(tmin, xinit, params)` : unused argument

The `pfun`function is probably missing one of the arguments that go into the function (either `t`, `X`or `params`). As an example:

```{r error3, eval = FALSE}
#THIS IS INCORRECT (MISSING T):
pfun       <- function(X,params){ cbind(params[1] *(1 + sin(t))* X[,1], 
                          (params[2] + (params[1]-params[2])*X[,1]/params[3])*X[,1] ) }

#THIS IS CORRECT:
pfun       <- function(t,X,params){ cbind(params[1] *(1 + sin(t))* X[,1], 
                          (params[2] + (params[1]-params[2])*X[,1]/params[3])*X[,1] ) }
```

### `nsim` is not a strictly positive integer. Defaulting to closest positive integer

If the number of simulations, `nsim` is smaller than 2, or is not an integer, the program automatically chooses a new value for  `nsim` looking for the closest integer. 

###  `kthsave` is not a strictly positive integer. Defaulting to closest positive integer

If the number indicating after how many iterations to save, `kthsave` is smaller than 2, or is not an integer, the program automatically chooses a new value for  `kthsave` looking for the closest positive integer. 

### `tmin >= tmax`

The time at which the simulation starts `tmin` is bigger or equal than the time at which the simulation ends `tmax`



## What is missing?

The project still needs a lot of testing and debugging. Furthermore, we are developing an automatic [tau-leaping algorithm](http://www.cs.ucsb.edu/~cse/Files/adaptivetau06.pdf)  to compliment our package. 

Please feel free to contribute to the project.  

## Contributor Code of Conduct

As contributors and maintainers of this project, we pledge to respect all people who contribute through reporting issues, posting feature requests, updating documentation, submitting pull requests or patches, and other activities.

We are committed to making participation in this project a harassment-free experience for everyone, regardless of level of experience, gender, gender identity and expression, sexual orientation, disability, personal appearance, body size, race, ethnicity, age, or religion.

Examples of unacceptable behavior by participants include the use of sexual language or imagery, derogatory comments or personal attacks, trolling, public or private harassment, insults, or other unprofessional conduct.

Project maintainers have the right and responsibility to remove, edit, or reject comments, commits, code, wiki edits, issues, and other contributions that are not aligned to this Code of Conduct. Project maintainers who do not follow the Code of Conduct may be removed from the project team.

Instances of abusive, harassing, or otherwise unacceptable behavior may be reported by opening an issue or contacting one or more of the project maintainers.

This Code of Conduct is adapted from the Contributor Covenant, version 1.0.0, available from [http://contributor-covenant.org/version/1/0/0/](http://contributor-covenant.org/version/1/0/0/)

## Licence

This package is free and open source software, licensed under [GPL-3](https://www.gnu.org/licenses/gpl-3.0.html). 

If you use this package please don't forget to cite it.

## Authors

* Rodrigo Zepeda-Tello [rodrigo.zepeda@insp.mx](mailto:rodrigo.zepeda@insp.mx)
* Dalia Camacho-García-Formentí [daliaf172@gmail.com](mailto:daliaf172@gmail.com)