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