benchmark_integrate_gamma.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Benchmarking integrators on statistical distributions from 3 langs. Matlab has a huge overhead"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "DIR = 'gammacdf'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.3678794411714423, 4.08428225874771e-15)\n",
      "CPU times: user 22.9 s, sys: 0 ns, total: 22.9 s\n",
      "Wall time: 22.9 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "#### Python\n",
    "import scipy.stats as spstat\n",
    "import scipy.integrate as spint\n",
    "f = spstat.gamma(a=1).cdf\n",
    "for i in range(int(1E4)):\n",
    "    b= spint.quad(f,0,1)\n",
    "print b\n",
    "\n",
    "# (0.3678794411714423, 4.08428225874771e-15)\n",
    "# CPU times: user 22.9 s, sys: 0 ns, total: 22.9 s\n",
    "# Wall time: 22.9 s\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.3678794411714423, 4.08428225874771e-15)\n",
      "CPU times: user 752 ms, sys: 24 ms, total: 776 ms\n",
      "Wall time: 743 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "import scipy.stats as spstat\n",
    "import scipy.integrate as spint\n",
    "import scipy.special as spspec\n",
    "\n",
    "f = lambda x: spstat.gamma._cdf(x,1)\n",
    "b= spint.quad(f,0,1)\n",
    "for i in range(int(1E4)):\n",
    "    b= spint.quad(f,0,1)\n",
    "print b\n",
    "\n",
    "# (0.3678794411714423, 4.08428225874771e-15)\n",
    "# CPU times: user 792 ms, sys: 16 ms, total: 808 ms\n",
    "# Wall time: 767 ms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  0.927352 seconds (1.34 M allocations: 35.109 MiB, 1.57% gc time)\n",
      "  0.296866 seconds (839.09 k allocations: 23.614 MiB, 1.31% gc time)\n",
      "(0.3678794411714423, 5.551115123125783e-17)\n",
      "\n",
      "CPU times: user 108 ms, sys: 16 ms, total: 124 ms\n",
      "Wall time: 5.82 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "### Julia\n",
    "lines = ''''\n",
    "using Distributions\n",
    "using Cubature\n",
    "D = Gamma(1)\n",
    "f = x ->begin cdf(D,x);end\n",
    "reps = 1E4\n",
    "@time for j in 1:reps \n",
    "    b = quadgk(f,0,1)\n",
    "    end\n",
    "@time for j in 1:reps \n",
    "    b = hquadrature(f,0,1)\n",
    "    end\n",
    "b = quadgk(f,0,1)\n",
    "println(b)' '''\n",
    "!echo {lines} > {DIR}/test.jl\n",
    "!julia {DIR}/test.jl\n",
    "print "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 152,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[?1h\u001b=\n",
      "                            < M A T L A B (R) >\n",
      "                  Copyright 1984-2016 The MathWorks, Inc.\n",
      "                   R2016b (9.1.0.441655) 64-bit (glnxa64)\n",
      "                             September 7, 2016\n",
      "\n",
      " \n",
      "To get started, type one of these: helpwin, helpdesk, or demo.\n",
      "For product information, visit www.mathworks.com.\n",
      " \n",
      "    0.3679\n",
      "\n",
      "test elapsed: 5.680399 s \u001b[?1l\u001b>\n",
      "\n",
      "CPU times: user 272 ms, sys: 36 ms, total: 308 ms\n",
      "Wall time: 13.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "### Matlab\n",
    "lines = r''' \"start = tic;\n",
    "f = @(x)gamcdf(x,1,1);\n",
    "for i = 1:1E4\n",
    "    b = integral(f,0,1);\n",
    "end\n",
    "disp(b)\n",
    "fprintf('%s elapsed: %f s ', mfilename,toc(start)) \" ''' \n",
    "!echo {lines} >{DIR}/test.m\n",
    "!matlab -nodisplay -r \"run('gammacdf/test.m');exit;\"\n",
    "print '\\n'\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 184,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   user  system elapsed \n",
      "  0.316   0.000   0.318 \n",
      "0.3678794 with absolute error < 4.1e-15\n",
      "\n",
      "CPU times: user 4 ms, sys: 20 ms, total: 24 ms\n",
      "Wall time: 634 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "### R\n",
    "lines = ''''\n",
    "main <- function(){\n",
    "    f = function(x){pgamma(x,1)}\n",
    "    reps = 1E4\n",
    "    for (j in 1:reps){\n",
    "        b = integrate(f,0,1)\n",
    "    }\n",
    "    b\n",
    "}\n",
    "system.time({b <- main()})\n",
    "print(b)\n",
    "' '''\n",
    "!echo {lines} > {DIR}/test.R\n",
    "!Rscript {DIR}/test.R\n",
    "print "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}