1  Benchmark between Python and Julia
1.1  The Romberg method
1.2  The Romberg method, naive recursive version in Python
1.3  The Romberg method, dynamic programming version in Python
1.4  The Romberg method, better dynamic programming version in Python
1.5  First benchmark
1.6  Using Pypy for speedup
1.7  Numba version for Python
1.8  Naive Julia version
1.9  Benchmark between Python, Pypy and Julia
1.10  Second benchmark
1.11  Conclusion
1.11.1  Remark
" ] }, { "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 = (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.5 * r[i - 1]) + 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\u001b[0m in \u001b[0;36m\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mromberg_numba\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m8\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;31m# Almost the exact value.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/dispatcher.py\u001b[0m in \u001b[0;36m_compile_for_args\u001b[0;34m(self, *args, **kws)\u001b[0m\n\u001b[1;32m 284\u001b[0m \u001b[0margtypes\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtypeof_pyval\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 285\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 286\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtuple\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margtypes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 287\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0merrors\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mTypingError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 288\u001b[0m \u001b[0;31m# Intercept typing error that may be due to an argument\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/dispatcher.py\u001b[0m in \u001b[0;36mcompile\u001b[0;34m(self, sig)\u001b[0m\n\u001b[1;32m 530\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 531\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cache_misses\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0msig\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 532\u001b[0;31m \u001b[0mcres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_compiler\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreturn_type\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 533\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0madd_overload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcres\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 534\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_cache\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msave_overload\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msig\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcres\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/dispatcher.py\u001b[0m in \u001b[0;36mcompile\u001b[0;34m(self, args, return_type)\u001b[0m\n\u001b[1;32m 79\u001b[0m \u001b[0mimpl\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 80\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mreturn_type\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mreturn_type\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 81\u001b[0;31m flags=flags, locals=self.locals)\n\u001b[0m\u001b[1;32m 82\u001b[0m \u001b[0;31m# Check typing error if object mode is used\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 83\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mcres\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtyping_error\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0menable_pyobject\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/compiler.py\u001b[0m in \u001b[0;36mcompile_extra\u001b[0;34m(typingctx, targetctx, func, args, return_type, flags, locals, library)\u001b[0m\n\u001b[1;32m 691\u001b[0m pipeline = Pipeline(typingctx, targetctx, library,\n\u001b[1;32m 692\u001b[0m args, return_type, flags, locals)\n\u001b[0;32m--> 693\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mpipeline\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompile_extra\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 694\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 695\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/compiler.py\u001b[0m in \u001b[0;36mcompile_extra\u001b[0;34m(self, func)\u001b[0m\n\u001b[1;32m 348\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlifted\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 349\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlifted_from\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 350\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_compile_bytecode\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 351\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 352\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mcompile_ir\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunc_ir\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlifted\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mlifted_from\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/compiler.py\u001b[0m in \u001b[0;36m_compile_bytecode\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 656\u001b[0m \"\"\"\n\u001b[1;32m 657\u001b[0m \u001b[0;32massert\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfunc_ir\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 658\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_compile_core\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 659\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 660\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_compile_ir\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/compiler.py\u001b[0m in \u001b[0;36m_compile_core\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 643\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 644\u001b[0m \u001b[0mpm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfinalize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 645\u001b[0;31m \u001b[0mres\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mpm\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstatus\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 646\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mres\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 647\u001b[0m \u001b[0;31m# Early pipeline completion\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/compiler.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, status)\u001b[0m\n\u001b[1;32m 234\u001b[0m \u001b[0;31m# No more fallback pipelines?\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 235\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mis_final_pipeline\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 236\u001b[0;31m \u001b[0;32mraise\u001b[0m \u001b[0mpatched_exception\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 237\u001b[0m \u001b[0;31m# Go to next fallback pipeline\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 238\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/compiler.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self, status)\u001b[0m\n\u001b[1;32m 226\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 227\u001b[0m \u001b[0mevent\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mstage_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 228\u001b[0;31m \u001b[0mstage\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 229\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0m_EarlyPipelineCompletion\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 230\u001b[0m \u001b[0;32mreturn\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mresult\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/compiler.py\u001b[0m in \u001b[0;36mstage_analyze_bytecode\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 362\u001b[0m \u001b[0mAnalyze\u001b[0m \u001b[0mbytecode\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mtranslating\u001b[0m \u001b[0mto\u001b[0m \u001b[0mNumba\u001b[0m \u001b[0mIR\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 363\u001b[0m \"\"\"\n\u001b[0;32m--> 364\u001b[0;31m \u001b[0mfunc_ir\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mtranslate_stage\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfunc_id\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 365\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_set_and_check_ir\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc_ir\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 366\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/compiler.py\u001b[0m in \u001b[0;36mtranslate_stage\u001b[0;34m(func_id, bytecode)\u001b[0m\n\u001b[1;32m 754\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mtranslate_stage\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc_id\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbytecode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 755\u001b[0m \u001b[0minterp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minterpreter\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mInterpreter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfunc_id\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 756\u001b[0;31m \u001b[0;32mreturn\u001b[0m \u001b[0minterp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0minterpret\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbytecode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 757\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 758\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/interpreter.py\u001b[0m in \u001b[0;36minterpret\u001b[0;34m(self, bytecode)\u001b[0m\n\u001b[1;32m 89\u001b[0m \u001b[0;31m# Control flow analysis\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 90\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcfa\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcontrolflow\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mControlFlowAnalysis\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mbytecode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 91\u001b[0;31m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcfa\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mrun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 92\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mconfig\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mDUMP_CFG\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 93\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcfa\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdump\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;32m/usr/local/lib/python3.5/dist-packages/numba/controlflow.py\u001b[0m in \u001b[0;36mrun\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m 513\u001b[0m \u001b[0mfn\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0minst\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 514\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 515\u001b[0;31m \u001b[0;32massert\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0minst\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_jump\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0minst\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 516\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 517\u001b[0m \u001b[0;31m# Close all blocks\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", "\u001b[0;31mAssertionError\u001b[0m: Failed at object (analyzing bytecode)\nSETUP_EXCEPT(arg=82, lineno=17)" ] } ], "source": [ "romberg_numba(f, a, b, n=8) # Almost the exact value." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> It fails! Almost as always when trying Numba, it fails cryptically, too bad. I don't want to spend time debugging this." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "## Naive Julia version" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "> Thanks to [this page](https://learnxinyminutes.com/docs/julia/) for a nice and short introduction to Julia." ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "code_folding": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "f (generic function with 1 method)\n", "1993\n", "2017\n", "romberg_julia (generic function with 1 method)\n", "404122.6272072803\n", "372300.32065714896\n", "373621.99733807985\n", "373650.47843482986\n", "409937.6105242601\n", "409937.57869815256\n", "409937.5786979804\n" ] } ], "source": [ "%%script julia\n", "\n", "function f(x)\n", " (12*x + 1) / (1 + cos(x)^2)\n", "end\n", "\n", "a = 1993\n", "b = 2017\n", "\n", "function romberg_julia(f, xmin, xmax; n=8)\n", " m = n\n", " # First value:\n", " r = Dict()\n", " r[(0, 0)] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.\n", "\n", " # One side of the triangle:\n", " for i in 1 : n\n", " h_i = (xmax - xmin) / (2^i)\n", " sum_f_x = 0\n", " for k in 1 : 2^(i - 1)\n", " sum_f_x += f(xmin + ((2 * k - 1) * h_i))\n", " end\n", " r[(i, 0)] = (r[(i - 1, 0)] / 2.) + (h_i * sum_f_x)\n", " end\n", "\n", " # All the other values:\n", " for j in 1 : m\n", " for i in j : n\n", " r[(i, j)] = (((4^j) * r[(i, j - 1)]) - r[(i - 1, j - 1)]) / (4^j - 1.)\n", " end\n", " end\n", "\n", " r[(n, m)]\n", "end\n", "\n", "\n", "println(romberg_julia(f, a, b, n=0)) # really not accurate!\n", "println(romberg_julia(f, a, b, n=1)) # alreay pretty good!\n", "println(romberg_julia(f, a, b, n=2))\n", "println(romberg_julia(f, a, b, n=3))\n", "println(romberg_julia(f, a, b, n=8)) # Almost the exact value.\n", "println(romberg_julia(f, a, b, n=10)) # Almost the exact value.\n", "println(romberg_julia(f, a, b, n=12)) # Almost the exact value.\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It seems to work well, like the Python implementation. We get the same numerical result:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(409937.57869797916, 8.482168070998719e-05)" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" }, { "data": { "text/plain": [ "409937.5786979791" ] }, "execution_count": 20, "metadata": {}, "output_type": "execute_result" } ], "source": [ "f = lambda x: (12*x+1)/(1+math.cos(x)**2)\n", "a, b = 1993, 2017\n", "\n", "quad(f, a, b)\n", "romberg(f, a, b, n=12)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let try a less naive version using a fixed-size array instead of a dictionary. (as we did before for the Python version)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": { "code_folding": [] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "f (generic function with 1 method)\n", "1993\n", "2017\n", "romberg_julia_better (generic function with 1 method)\n", "404122.6272072803\n", "372300.32065714896\n", "373621.99733807985\n", "373650.47843482986\n", "409937.6105242601\n", "409937.57869815256\n", "409937.5786979804\n" ] } ], "source": [ "%%script julia\n", "\n", "function f(x)\n", " (12*x + 1) / (1 + cos(x)^2)\n", "end\n", "\n", "a = 1993\n", "b = 2017\n", "\n", "function romberg_julia_better(f, xmin, xmax; n=8)\n", " m = n\n", " # First value:\n", " r = zeros((n+1, m+1)) # https://docs.julialang.org/en/latest/stdlib/arrays/#Base.zeros\n", " r[1, 1] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.\n", "\n", " # One side of the triangle:\n", " for i in 1 : n\n", " h_i = (xmax - xmin) / (2^i)\n", " sum_f_x = 0\n", " for k in 1 : 2^(i - 1)\n", " sum_f_x += f(xmin + ((2 * k - 1) * h_i))\n", " end\n", " r[i + 1, 1] = (r[i, 1] / 2.) + (h_i * sum_f_x)\n", " end\n", "\n", " # All the other values:\n", " for j in 1 : m\n", " for i in j : n\n", " r[i + 1, j + 1] = (((4.^j) * r[i + 1, j]) - r[i, j]) / (4.^j - 1.)\n", " end\n", " end\n", "\n", " r[n + 1, m + 1]\n", "end\n", "\n", "\n", "println(romberg_julia_better(f, a, b, n=0)) # really not accurate!\n", "println(romberg_julia_better(f, a, b, n=1)) # alreay pretty good!\n", "println(romberg_julia_better(f, a, b, n=2))\n", "println(romberg_julia_better(f, a, b, n=3))\n", "println(romberg_julia_better(f, a, b, n=8)) # Almost the exact value.\n", "println(romberg_julia_better(f, a, b, n=10)) # Almost the exact value.\n", "println(romberg_julia_better(f, a, b, n=12)) # Almost the exact value." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "## Benchmark between Python, Pypy and Julia" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First with Python:" ] }, { "cell_type": "code", "execution_count": 22, "metadata": { "code_folding": [ 9 ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 25.6 s, sys: 0 ns, total: 25.6 s\n", "Wall time: 25.6 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": 23, "metadata": { "code_folding": [ 9 ], "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CPU times: user 4 ms, sys: 4 ms, total: 8 ms\n", "Wall time: 7.12 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 = (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.5 * r[i - 1]) + 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": [ "And finally with Julia:" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "code_folding": [ 7 ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "f (generic function with 1 method)\n", "romberg_julia (generic function with 1 method)\n", "CPU times: user 4 ms, sys: 4 ms, total: 8 ms\n", "Wall time: 8.1 s\n" ] } ], "source": [ "%%time\n", "%%script julia\n", "\n", "function f(x)\n", " (12*x + 1) / (1 + cos(x)^2)\n", "end\n", "\n", "function romberg_julia(f, xmin, xmax; n=8)\n", " m = n\n", " # First value:\n", " r = Dict()\n", " r[(0, 0)] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.\n", "\n", " # One side of the triangle:\n", " for i in 1 : n\n", " h_i = (xmax - xmin) / (2^i)\n", " sum_f_x = 0\n", " for k in 1 : 2^(i - 1)\n", " sum_f_x += f(xmin + ((2 * k - 1) * h_i))\n", " end\n", " r[(i, 0)] = (r[(i - 1, 0)] / 2.) + (h_i * sum_f_x)\n", " end\n", "\n", " # All the other values:\n", " for j in 1 : m\n", " for i in j : n\n", " r[(i, j)] = (((4^j) * r[(i, j - 1)]) - r[(i - 1, j - 1)]) / (4^j - 1.)\n", " end\n", " end\n", "\n", " r[(n, m)]\n", "end\n", "\n", "for _ in 1:100000\n", " a = rand(-2000:2000)\n", " b = a + rand(0:100)\n", " romberg_julia(f, a, b)\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "On this first test, it doesn't look faster than Pypy...\n", "But what if we use the improved version, with an array instead of dictionary?" ] }, { "cell_type": "code", "execution_count": 25, "metadata": { "code_folding": [ 7 ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "f (generic function with 1 method)\n", "romberg_julia_better (generic function with 1 method)\n", "CPU times: user 4 ms, sys: 4 ms, total: 8 ms\n", "Wall time: 2.7 s\n" ] } ], "source": [ "%%time\n", "%%script julia\n", "\n", "function f(x)\n", " (12*x + 1) / (1 + cos(x)^2)\n", "end\n", "\n", "function romberg_julia_better(f, xmin, xmax; n=8)\n", " m = n\n", " # First value:\n", " r = zeros((n+1, m+1)) # https://docs.julialang.org/en/latest/stdlib/arrays/#Base.zeros\n", " r[1, 1] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.\n", "\n", " # One side of the triangle:\n", " for i in 1 : n\n", " h_i = (xmax - xmin) / (2^i)\n", " sum_f_x = 0\n", " for k in 1 : 2^(i - 1)\n", " sum_f_x += f(xmin + ((2 * k - 1) * h_i))\n", " end\n", " r[i + 1, 1] = (r[i, 1] / 2.) + (h_i * sum_f_x)\n", " end\n", "\n", " # All the other values:\n", " for j in 1 : m\n", " for i in j : n\n", " r[i + 1, j + 1] = (((4.^j) * r[i + 1, j]) - r[i, j]) / (4.^j - 1.)\n", " end\n", " end\n", "\n", " r[n + 1, m + 1]\n", "end\n", "\n", "for _ in 1:100000\n", " a = rand(-2000:2000)\n", " b = a + rand(0:100)\n", " romberg_julia_better(f, a, b)\n", "end" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Oh, this time it finally seems faster. Really faster? Yes, about 3 to 4 time faster than Pypy.\n", "\n", "Remark also that this last cells compared by using the magic %%pypy and %%script julia, so they both need a warm-up time (opening the pipe, the sub-process, initializing the JIT compiler etc).\n", "But it is fair to compare Pypy to Julia this way." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "## Second benchmark" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let try the same numerical algorithm but with a different integrand function.\n", "\n", "$$\\int_{0}^{1} \\exp(-x^2) \\mathrm{d}x \\approx 0.842700792949715$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First with Python:" ] }, { "cell_type": "code", "execution_count": 26, "metadata": { "code_folding": [ 9 ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.84270079295\n", "CPU times: user 20.7 s, sys: 8 ms, total: 20.8 s\n", "Wall time: 20.8 s\n" ] } ], "source": [ "%%time\n", "\n", "import numpy as np\n", "import math\n", "import random\n", "\n", "f = lambda x: (2.0 / math.sqrt(math.pi)) * math.exp(-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 = 0\n", " b = 1\n", " romberg(f, a, b)\n", "print(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": 27, "metadata": { "code_folding": [ 9 ], "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.84270079295\n", "CPU times: user 8 ms, sys: 0 ns, total: 8 ms\n", "Wall time: 6.35 s\n" ] } ], "source": [ "%%time\n", "%%pypy\n", "\n", "import math\n", "import random\n", "\n", "f = lambda x: (2.0 / math.sqrt(math.pi)) * math.exp(-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 = (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.5 * r[i - 1]) + 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 = 0\n", " b = 1\n", " romberg_pypy(f, a, b)\n", "print(romberg_pypy(f, a, b))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "And finally with Julia:" ] }, { "cell_type": "code", "execution_count": 28, "metadata": { "code_folding": [ 7 ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "f (generic function with 1 method)\n", "romberg_julia (generic function with 1 method)\n", "0.8427007929497148\n", "CPU times: user 4 ms, sys: 4 ms, total: 8 ms\n", "Wall time: 7.19 s\n" ] } ], "source": [ "%%time\n", "%%script julia\n", "\n", "function f(x)\n", " (2.0 / sqrt(pi)) * exp(-x^2)\n", "end\n", "\n", "function romberg_julia(f, xmin, xmax; n=8)\n", " m = n\n", " # First value:\n", " r = Dict()\n", " r[(0, 0)] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.\n", "\n", " # One side of the triangle:\n", " for i in 1 : n\n", " h_i = (xmax - xmin) / (2^i)\n", " sum_f_x = 0\n", " for k in 1 : 2^(i - 1)\n", " sum_f_x += f(xmin + ((2 * k - 1) * h_i))\n", " end\n", " r[(i, 0)] = (r[(i - 1, 0)] / 2.) + (h_i * sum_f_x)\n", " end\n", "\n", " # All the other values:\n", " for j in 1 : m\n", " for i in j : n\n", " r[(i, j)] = (((4^j) * r[(i, j - 1)]) - r[(i - 1, j - 1)]) / (4^j - 1.)\n", " end\n", " end\n", "\n", " r[(n, m)]\n", "end\n", "\n", "for _ in 1:100000\n", " a = 0\n", " b = 1\n", " romberg_julia(f, a, b)\n", "end\n", "println(romberg_julia(f, 0, 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Still not faster than Pypy... So what is the goal of Julia?" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "code_folding": [ 8 ] }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "f (generic function with 1 method)\n", "romberg_julia_better (generic function with 1 method)\n", "0.8427007929497148\n", "CPU times: user 8 ms, sys: 0 ns, total: 8 ms\n", "Wall time: 2.42 s\n" ] } ], "source": [ "%%time\n", "%%script julia\n", "\n", "function f(x)\n", " (2.0 / sqrt(pi)) * exp(-x^2)\n", "end\n", "\n", "\n", "function romberg_julia_better(f, xmin, xmax; n=8)\n", " m = n\n", " # First value:\n", " r = zeros((n+1, m+1)) # https://docs.julialang.org/en/latest/stdlib/arrays/#Base.zeros\n", " r[1, 1] = (xmax - xmin) * (f(xmax) + f(xmin)) / 2.\n", "\n", " # One side of the triangle:\n", " for i in 1 : n\n", " h_i = (xmax - xmin) / (2^i)\n", " sum_f_x = 0\n", " for k in 1 : 2^(i - 1)\n", " sum_f_x += f(xmin + ((2 * k - 1) * h_i))\n", " end\n", " r[i + 1, 1] = (r[i, 1] / 2.) + (h_i * sum_f_x)\n", " end\n", "\n", " # All the other values:\n", " for j in 1 : m\n", " for i in j : n\n", " r[i + 1, j + 1] = (((4.^j) * r[i + 1, j]) - r[i, j]) / (4.^j - 1.)\n", " end\n", " end\n", "\n", " r[n + 1, m + 1]\n", "end\n", "\n", "for _ in 1:100000\n", " a = 0\n", " b = 1\n", " romberg_julia_better(f, a, b)\n", "end\n", "println(romberg_julia_better(f, 0, 1))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This is also faster than Pypy, but not that much..." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "----\n", "## Conclusion\n", "\n", "$\\implies$\n", "On this (baby) example of a real world numerical algorithms, tested on thousands of random inputs or on thousands time the same input, the speed-up is in favor of Julia, but it doesn't seem impressive enough to make me want to use it (at least for now).\n", "\n", "If I have to use 1-based indexing and a slightly different language, just to maybe gain a speed-up of 2 to 3 (compared to Pypy) or even a 10x speed-up compare to naive Python, why bother?\n", "\n", "### Remark\n", "Of course, this was a *baby* benchmark, on a small algorithm, and probably wrongly implemented in both Python and Julia.\n", "\n", "But still, I am surprised to see that the naive Julia version was *slower* than the naive Python version executed with Pypy...\n", "For the less naive version (without dictionary), the Julia version was about *2 to 3* times faster than the Python version with Pypy." ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.5.3" }, "toc": { "colors": { "hover_highlight": "#DAA520", "running_highlight": "#FF0000", "selected_highlight": "#FFD700" }, "moveMenuLeft": true, "nav_menu": { "height": "258px", "width": "251px" }, "navigate_menu": true, "number_sections": true, "sideBar": true, "threshold": 4, "toc_cell": true, "toc_position": { "height": "506px", "left": "0px", "right": "1068px", "top": "116px", "width": "212px" }, "toc_section_display": "block", "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }