{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Benchmark between Python and Julia\n", "\n", "This small [Jupyter notebook](http://jupyter.org/) shows a simple benchmark comparing various implementations in Python and one in Julia of a specific numerical algorithm, the [Romberg integration method](https://en.wikipedia.org/wiki/Romberg%27s_method).\n", "\n", "For Python:\n", "\n", "- a recursive implementation,\n", "- a dynamic programming implementation,\n", "- also using Pypy instead,\n", "- (maybe a Numba version of the dynamic programming version)\n", " + (maybe a Cython version too)\n", "\n", "For Julia:\n", "\n", "- a dynamic programming implementation will be enough." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "## The Romberg method\n", "\n", "> For mathematical explanations, see [the Wikipedia page](https://en.wikipedia.org/wiki/Romberg%27s_method)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will use [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html) function to compare the result of our manual implementations." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from scipy.integrate import quad" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let try it with this function $f(x)$ on $[a,b]=[1993,2015]$:\n", "\n", "$$ f(x) := \\frac{12x+1}{1+\\cos(x)^2} $$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import math\n", "\n", "f = lambda x: (12*x+1)/(1+math.cos(x)**2)\n", "a, b = 1993, 2017" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(409937.57869797916, 8.482168070998719e-05)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "quad(f, a, b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first value is the numerical value of the integral $\\int_{a}^{b} f(x) \\mathrm{d}x$ and the second value is an estimate of the numerical error.\n", "\n", "$0.4\\%$ is not much, alright." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "## The Romberg method, naive recursive version in Python\n", "\n", "See https://mec-cs101-integrals.readthedocs.io/en/latest/_modules/integrals.html#romberg_rec for the code and https://mec-cs101-integrals.readthedocs.io/en/latest/integrals.html#integrals.romberg_rec for the doc" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def romberg_rec(f, xmin, xmax, n=8, m=None):\n", " if m is None: # not m was considering 0 as None\n", " m = n\n", " assert n >= m\n", " if n == 0 and m == 0:\n", " return ((xmax - xmin) / 2.0) * (f(xmin) + f(xmax))\n", " elif m == 0:\n", " h = (xmax - xmin) / float(2**n)\n", " N = (2**(n - 1)) + 1\n", " term = math.fsum(f(xmin + ((2 * k) - 1) * h) for k in range(1, N))\n", " return (term * h) + (0.5) * romberg_rec(f, xmin, xmax, n - 1, 0)\n", " else:\n", " return (1.0 / ((4**m) - 1)) * ((4**m) * romberg_rec(f, xmin, xmax, n, m - 1) - romberg_rec(f, xmin, xmax, n - 1, m - 1))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "404122.6272072803" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "372300.32065714896" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "373621.9973380798" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "373650.4784348298" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "409937.6105242601" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "409937.57869815244" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "409937.57869797904" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "romberg_rec(f, a, b, n=0) # really not accurate!\n", "romberg_rec(f, a, b, n=1) # alreay pretty good!\n", "romberg_rec(f, a, b, n=2)\n", "romberg_rec(f, a, b, n=3)\n", "romberg_rec(f, a, b, n=8) # Almost the exact value.\n", "romberg_rec(f, a, b, n=10) # Almost the exact value.\n", "romberg_rec(f, a, b, n=12) # Almost the exact value." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It converges quite quickly to the \"true\" value as given by `scipy.integrate.quad`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "## The Romberg method, dynamic programming version in Python\n", "\n", "See https://mec-cs101-integrals.readthedocs.io/en/latest/_modules/integrals.html#romberg for the code and https://mec-cs101-integrals.readthedocs.io/en/latest/integrals.html#integrals.romberg for the doc." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It is not hard to make this function non-recursive, by storing the intermediate results." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def romberg(f, xmin, xmax, n=8, m=None):\n", " assert xmin <= xmax\n", " if m is None:\n", " m = n\n", " assert n >= m >= 0\n", " # First value:\n", " r = {(0, 0): 0.5 * (xmax - xmin) * (f(xmax) + f(xmin))}\n", "\n", " # One side of the triangle:\n", " for i in range(1, n + 1):\n", " h_i = (xmax - xmin) / float(2**i)\n", " xsamples = [xmin + ((2 * k - 1) * h_i) for k in range(1, 1 + 2**(i - 1))]\n", " r[(i, 0)] = (0.5 * r[(i - 1, 0)]) + h_i * math.fsum(f(x) for x in xsamples)\n", "\n", " # All the other values:\n", " for j in range(1, m + 1):\n", " for i in range(j, n + 1):\n", " try:\n", " r[(i, j)] = (((4**j) * r[(i, j - 1)]) - r[(i - 1, j - 1)]) / float((4**j) - 1)\n", " except:\n", " raise ValueError(\"romberg() with n = {}, m = {} and i = {}, j = {} was an error.\".format(n, m, i, j))\n", "\n", " return r[(n, m)]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "404122.6272072803" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "372300.32065714896" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "373621.99733807985" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "373650.47843482986" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "409937.61052426015" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "409937.5786981525" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "409937.5786979791" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "romberg(f, a, b, n=0) # really not accurate!\n", "romberg(f, a, b, n=1) # alreay pretty good!\n", "romberg(f, a, b, n=2)\n", "romberg(f, a, b, n=3)\n", "romberg(f, a, b, n=8) # Almost the exact value.\n", "romberg(f, a, b, n=10) # Almost the exact value.\n", "romberg(f, a, b, n=12) # Almost the exact value." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It converges quite quickly to the \"true\" value as given by `scipy.integrate.quad`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "## The Romberg method, better dynamic programming version in Python\n", "\n", "Instead of using a dictionary, which gets filled up dynamically (and so, slowly), let us use a numpy arrays, as we already know the size of the array we need ($n+1 \\times m+1$).\n", "\n", "Note that only half of the array is used, so we could try to use [sparse matrices](https://docs.scipy.org/doc/scipy/reference/sparse.html) maybe, for triangular matrices? From what I know, it is not worth it, both in term of memory information (if the sparsity measure is $\\simeq 1/2$, you don't win anything from [LIL](https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.lil_matrix.html#scipy.sparse.lil_matrix) or other sparse matrices representation...\n", "We could use [`numpy.tri`](https://docs.scipy.org/doc/numpy/reference/generated/numpy.tri.html) but this uses a dense array so nope." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true }, "outputs": [], "source": [ "import numpy as np\n", "\n", "def romberg_better(f, xmin, xmax, n=8, m=None):\n", " assert xmin <= xmax\n", " if m is None:\n", " m = n\n", " assert n >= m >= 0\n", " # First value:\n", " r = np.zeros((n+1, m+1))\n", " r[0, 0] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.\n", "\n", " # One side of the triangle:\n", " for i in range(1, n + 1):\n", " h_i = (xmax - xmin) / 2.**i\n", " r[i, 0] = (0.5 * r[i - 1, 0]) + h_i * math.fsum(\n", " f(xmin + ((2 * k - 1) * h_i))\n", " for k in range(1, 1 + 2**(i - 1))\n", " )\n", "\n", " # All the other values:\n", " for j in range(1, m + 1):\n", " for i in range(j, n + 1):\n", " r[i, j] = (((4.**j) * r[i, j - 1]) - r[i - 1, j - 1]) / ((4.**j) - 1.)\n", "\n", " return r[n, m]" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "404122.62720728031" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "372300.32065714896" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "373621.99733807985" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "373650.47843482986" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "409937.61052426015" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "409937.5786981525" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "409937.5786979791" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "romberg_better(f, a, b, n=0) # really not accurate!\n", "romberg_better(f, a, b, n=1) # alreay pretty good!\n", "romberg_better(f, a, b, n=2)\n", "romberg_better(f, a, b, n=3)\n", "romberg_better(f, a, b, n=8) # Almost the exact value.\n", "romberg_better(f, a, b, n=10) # Almost the exact value.\n", "romberg_better(f, a, b, n=12) # Almost the exact value." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It converges quite quickly to the \"true\" value as given by `scipy.integrate.quad`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "## First benchmark" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "390 µs ± 10.7 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n" ] } ], "source": [ "%timeit quad(f, a, b)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "52.8 ms ± 1.21 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n" ] } ], "source": [ "%timeit romberg_rec(f, a, b, n=10)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "819 µs ± 19.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n" ] } ], "source": [ "%timeit romberg(f, a, b, n=10)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "844 µs ± 19 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n" ] } ], "source": [ "%timeit romberg_better(f, a, b, n=10)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We already see that the recursive version is *much* slower than the dynamic programming one!\n", "\n", "But there is not much difference between the one using dictionary (`romberg()`) and the one using a numpy array of a known size (`romberg_better()`)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "## Using Pypy for speedup" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "code_folding": [ 9 ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 24.1 s, sys: 0 ns, total: 24.1 s\n", "Wall time: 24.1 s\n" ] } ], "source": [ "%%time\n", "\n", "import numpy as np\n", "import math\n", "import random\n", "\n", "f = lambda x: (12*x+1)/(1+math.cos(x)**2)\n", "\n", "# Same code\n", "def romberg(f, xmin, xmax, n=8, m=None):\n", " assert xmin <= xmax\n", " if m is None:\n", " m = n\n", " assert n >= m >= 0\n", " # First value:\n", " r = np.zeros((n+1, m+1))\n", " r[0, 0] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.\n", "\n", " # One side of the triangle:\n", " for i in range(1, n + 1):\n", " h_i = (xmax - xmin) / 2.**i\n", " r[i, 0] = (0.5 * r[i - 1, 0]) + h_i * math.fsum(\n", " f(xmin + ((2 * k - 1) * h_i))\n", " for k in range(1, 1 + 2**(i - 1))\n", " )\n", "\n", " # All the other values:\n", " for j in range(1, m + 1):\n", " for i in range(j, n + 1):\n", " r[i, j] = (((4.**j) * r[i, j - 1]) - r[i - 1, j - 1]) / ((4.**j) - 1.)\n", "\n", " return r[n, m]\n", "\n", "for _ in range(100000):\n", " a = random.randint(-2000, 2000)\n", " b = a + random.randint(0, 100)\n", " romberg(f, a, b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And now the same code executed by an external [Pypy](http://pypy.org) interpreter (Python 2.7.13 and PyPy 5.8.0 with GCC 5.4.0)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "code_folding": [ 9 ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 4 ms, sys: 0 ns, total: 4 ms\n", "Wall time: 7.09 s\n" ] } ], "source": [ "%%time\n", "%%pypy\n", "\n", "import math\n", "import random\n", "\n", "f = lambda x: (12*x+1)/(1+math.cos(x)**2)\n", "\n", "# Same code\n", "def romberg_pypy(f, xmin, xmax, n=8, m=None):\n", " assert xmin <= xmax\n", " if m is None:\n", " m = n\n", " assert n >= m >= 0\n", " # First value:\n", " r = [[0 for _ in range(n+1)] for _ in range(m+1)]\n", " r[0][0] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.\n", "\n", " # One side of the triangle:\n", " for i in range(1, n + 1):\n", " h_i = (xmax - xmin) / 2.**i\n", " r[i][0] = (0.5 * r[i - 1][0]) + h_i * math.fsum(\n", " f(xmin + ((2 * k - 1) * h_i))\n", " for k in range(1, 1 + 2**(i - 1))\n", " )\n", "\n", " # All the other values:\n", " for j in range(1, m + 1):\n", " for i in range(j, n + 1):\n", " r[i][j] = (((4.**j) * r[i][j - 1]) - r[i - 1][j - 1]) / ((4.**j) - 1.)\n", "\n", " return r[n][m]\n", "\n", "for _ in range(100000):\n", " a = random.randint(-2000, 2000)\n", " b = a + random.randint(0, 100)\n", " romberg_pypy(f, a, b)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> This version uses the improved memoization trick (no dictionary), but uses nested lists and not numpy arrays, I didn't bother to install numpy on my Pypy installation (even thought [it should be possible](https://bitbucket.org/pypy/numpy.git))." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "## Numba version for Python" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from numba import jit" ] }, { "cell_type": "code", "execution_count": 17, "metadata": { "collapsed": true }, "outputs": [], "source": [ "@jit\n", "def romberg_numba(f, xmin, xmax, n=8):\n", " assert xmin <= xmax\n", " m = n\n", " # First value:\n", " r = {(0, 0): 0.5 * (xmax - xmin) * (f(xmax) + f(xmin))}\n", "\n", " # One side of the triangle:\n", " for i in range(1, n + 1):\n", " h_i = (xmax - xmin) / float(2**i)\n", " xsamples = [xmin + ((2 * k - 1) * h_i) for k in range(1, 1 + 2**(i - 1))]\n", " r[(i, 0)] = (0.5 * r[(i - 1, 0)]) + h_i * math.fsum(f(x) for x in xsamples)\n", "\n", " # All the other values:\n", " for j in range(1, m + 1):\n", " for i in range(j, n + 1):\n", " try:\n", " r[(i, j)] = (((4**j) * r[(i, j - 1)]) - r[(i - 1, j - 1)]) / float((4**j) - 1)\n", " except:\n", " raise ValueError(\"romberg() with n = {}, m = {} and i = {}, j = {} was an error.\".format(n, m, i, j))\n", "\n", " return r[(n, m)]" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [ { "ename": "AssertionError", "evalue": "Failed at object (analyzing bytecode)\nSETUP_EXCEPT(arg=82, lineno=17)", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mAssertionError\u001b[0m Traceback (most recent call last)", "\u001b[0;32m