diff --git a/Benchmark_between_Python_and_Julia.ipynb b/Benchmark_between_Python_and_Julia.ipynb
new file mode 100644
index 0000000..f1f4fd2
--- /dev/null
+++ b/Benchmark_between_Python_and_Julia.ipynb
@@ -0,0 +1,1290 @@
+{
+ "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",
+ "## First benchmark"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "372 µs ± 11 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
+ ]
+ }
+ ],
+ "source": [
+ "%timeit quad(f, a, b)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "51.9 ms ± 803 µs 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": 10,
+ "metadata": {},
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "801 µs ± 8.32 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n"
+ ]
+ }
+ ],
+ "source": [
+ "%timeit romberg(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!"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "----\n",
+ "## Using Pypy for speedup"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 11,
+ "metadata": {
+ "code_folding": [
+ 7
+ ]
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "CPU times: user 27.8 s, sys: 4 ms, total: 27.8 s\n",
+ "Wall time: 27.8 s\n"
+ ]
+ }
+ ],
+ "source": [
+ "%%time\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(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)]\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": 12,
+ "metadata": {
+ "code_folding": [
+ 9
+ ]
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "CPU times: user 0 ns, sys: 8 ms, total: 8 ms\n",
+ "Wall time: 7.84 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, 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)]\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": [
+ "----\n",
+ "## Numba version for Python"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 13,
+ "metadata": {
+ "collapsed": true
+ },
+ "outputs": [],
+ "source": [
+ "from numba import jit"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 14,
+ "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": 15,
+ "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!"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "----\n",
+ "## Naive Julia version"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "> Thanks to [this page](https://learnxinyminutes.com/docs/julia/)."
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 16,
+ "metadata": {},
+ "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": 17,
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "(409937.57869797916, 8.482168070998719e-05)"
+ ]
+ },
+ "execution_count": 17,
+ "metadata": {},
+ "output_type": "execute_result"
+ },
+ {
+ "data": {
+ "text/plain": [
+ "409937.5786979791"
+ ]
+ },
+ "execution_count": 17,
+ "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": [
+ "----\n",
+ "## Benchmark between Python, Pypy and Julia"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "metadata": {},
+ "source": [
+ "First with Python:"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 18,
+ "metadata": {
+ "code_folding": [
+ 8
+ ]
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "CPU times: user 26.1 s, sys: 0 ns, total: 26.1 s\n",
+ "Wall time: 26.1 s\n"
+ ]
+ }
+ ],
+ "source": [
+ "%%time\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(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)]\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": 19,
+ "metadata": {
+ "code_folding": [
+ 9
+ ],
+ "scrolled": true
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "CPU times: user 0 ns, sys: 4 ms, total: 4 ms\n",
+ "Wall time: 7.06 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, 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)]\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": 20,
+ "metadata": {
+ "code_folding": [
+ 9
+ ]
+ },
+ "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.32 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..."
+ ]
+ },
+ {
+ "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": 21,
+ "metadata": {
+ "code_folding": [
+ 8
+ ]
+ },
+ "outputs": [
+ {
+ "name": "stdout",
+ "output_type": "stream",
+ "text": [
+ "0.842700792949715\n",
+ "CPU times: user 22.1 s, sys: 4 ms, total: 22.1 s\n",
+ "Wall time: 22.1 s\n"
+ ]
+ }
+ ],
+ "source": [
+ "%%time\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(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)]\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": 22,
+ "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.34 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, 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)]\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": 23,
+ "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: 8 ms, total: 12 ms\n",
+ "Wall time: 7.07 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": "markdown",
+ "metadata": {},
+ "source": [
+ "----\n",
+ "## Conclusion\n",
+ "\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..."
+ ]
+ }
+ ],
+ "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_section_display": "block",
+ "toc_window_display": true
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}