Benchmarks for SIR model in R/Rcpp

sir-r-benchmark.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "library(GillespieSSA)\n",
    "library(Rcpp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "sir <- function(beta, gamma, N, S0, I0, R0, tf) {\n",
    "    time <- 0\n",
    "    S <- S0\n",
    "    I <- I0\n",
    "    R <- R0\n",
    "    ta <- numeric(0)\n",
    "    Sa <- numeric(0)\n",
    "    Ia <- numeric(0)\n",
    "    Ra <- numeric(0)\n",
    "    while (time < tf) {\n",
    "        ta <- c(ta, time)\n",
    "        Sa <- c(Sa, S)\n",
    "        Ia <- c(Ia, I)\n",
    "        Ra <- c(Ra, R)\n",
    "        pf1 <- beta * S * I\n",
    "        pf2 <- gamma * I\n",
    "        pf <- pf1 + pf2\n",
    "        dt <- rexp(1, rate = pf)\n",
    "        time <- time + dt\n",
    "        if (time > tf) {\n",
    "            break\n",
    "        }\n",
    "        ru <- runif(1)\n",
    "        if (ru < (pf1/pf)) {\n",
    "            S <- S - 1\n",
    "            I <- I + 1\n",
    "        } else {\n",
    "            I <- I - 1\n",
    "            R <- R + 1\n",
    "        }\n",
    "        if (I == 0) {\n",
    "            break\n",
    "        }\n",
    "    }\n",
    "    results <- data.frame(time = ta, S = Sa, I = Ia, R = Ra)\n",
    "    return(results)\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "cppFunction('\n",
    "  List sirc(double beta, double gamma, double N, double S0, double I0, double R0, double tf){\n",
    "    double t = 0;\n",
    "    double S = S0;\n",
    "    double I = I0;\n",
    "    double R = R0;\n",
    "    std::vector<double> ta;\n",
    "    std::vector<double> Sa;\n",
    "    std::vector<double> Ia;\n",
    "    std::vector<double> Ra;\n",
    "    do{\n",
    "      ta.push_back(t);\n",
    "      Sa.push_back(S);\n",
    "      Ia.push_back(I);\n",
    "      Ra.push_back(R);\n",
    "      double pf1 = beta*S*I;\n",
    "      double pf2 = gamma*I;\n",
    "      double pf = pf1+pf2;\n",
    "      double dt = rexp(1,pf)[0];\n",
    "      t += dt;\n",
    "      double r = runif(1)[0];\n",
    "      if(r<pf1/pf){\n",
    "        S--;\n",
    "        I++;\n",
    "      }else{\n",
    "        I--;\n",
    "        R++;\n",
    "      }\n",
    "      if(I==0){break;}\n",
    "    } while (t<=tf && (I>0));\n",
    "  return List::create(_[\"time\"] = ta, _[\"S\"] = Sa, _[\"I\"] = Ia, _[\"R\"]=Ra);\n",
    "  }'\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "N <- 1000\n",
    "# R\n",
    "set.seed(4)\n",
    "sirsim <- sum(replicate(N, system.time(x <- sir(0.1/10000, 0.05, 10000, 9999, 1, 0, 1000))[\"elapsed\"]), trim = 0.05)\n",
    "# Rcpp\n",
    "set.seed(4)\n",
    "sircsim <- sum(replicate(N, system.time(x <- sirc(0.1/10000, 0.05, 10000, 9999, 1, 0, 1000))[\"elapsed\"]), trim = 0.05)\n",
    "# GillespieSSA\n",
    "sir.parms <- c(beta=0.1/10000,gamma=0.05)\n",
    "sir.x0 <- c(S=9999,I=1,R=0)\n",
    "sir.a <- c(\"beta*S*I\",\"gamma*I\")\n",
    "sir.nu <- matrix(c(-1,+1,0,0,-1,+1),nrow=3,ncol=2,byrow=FALSE)\n",
    "tf <- 1000.0\n",
    "sirssasim <- sum(replicate(N, system.time(x <- ssa(sir.x0,sir.a,sir.nu,sir.parms,tf=tf,simName=\"SIR\"))[\"elapsed\"]), trim = 0.05)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[1] \"R: 1087.945 ms\"\n",
      "[1] \"Rcpp: 1.30500000000079 ms\"\n",
      "[1] \"GillespieSSA: 894.246999999999 ms\"\n"
     ]
    }
   ],
   "source": [
    "print(paste(\"R:\",sirsim/N*1000.0,\"ms\"))\n",
    "print(paste(\"Rcpp:\",sircsim/N*1000,\"ms\"))\n",
    "print(paste(\"GillespieSSA:\",sirssasim/N*1000,\"ms\"))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "R",
   "language": "R",
   "name": "ir"
  },
  "language_info": {
   "codemirror_mode": "r",
   "file_extension": ".r",
   "mimetype": "text/x-r-source",
   "name": "R",
   "pygments_lexer": "r",
   "version": "3.3.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}