{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"\n",
" \n",
"
Is it a bird?
Is it a plane?
Accelerating Python with numba
\n",
" Juan Luis Cano Rodríguez <hello@juanlu.space>
2020-04-09 @ PyAmsterdam #StayAtHome
\n",
" \n",
" "
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Outline\n",
"\n",
"1. \"Python is slow\"\n",
"2. What is numba?\n",
"3. Some examples\n",
"4. Limitations and workarounds\n",
"5. Conclusions"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# Who is this guy?\n",
"\n",
"* **Aerospace Engineer** with a passion for orbits 🛰\n",
"* Former chair of the **Python España** non profit and former co-organizer of **PyCon Spain** 🐍\n",
"* **Mission Planning & Execution Engineer** at **Satellogic** 🌍\n",
"* Free Software advocate and Python enthusiast 🕮\n",
"* Hard Rock lover 🎸\n",
"\n",
"Follow me! https://github.com/astrojuanlu/\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# \"Python is slow\"\n",
"\n",
"## Data structures\n",
"\n",
"\n",
"\n",
"_(From https://jakevdp.github.io/blog/2014/05/09/why-python-is-slow/_)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Dynamic and interpreted (rather than static and compiled)\n",
"\n",
"\n",
"\n",
"_(From https://gist.github.com/Juanlu001/cf19b1c16caf618860fb_)"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Introspection\n",
"\n",
"\n",
"\n",
"_(From https://stackoverflow.com/q/42378726/554319)_"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## What to do?\n",
"\n",
"### Vectorization\n",
"\n",
"* Rewriting some code leveraging high level NumPy functions can make it way faster\n",
"* However, this works best for array manipulation - some other algorithms cannot easily be vectorized\n",
"* And even if you can, vectorized code can be impossible to read\n",
"\n",
""
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### Cython\n",
"\n",
"\n",
"\n",
"* Mature, widely used, effective, gradual - a great project!\n",
"* Some personal problems with it:\n",
" - I don't know any C, so it's more difficult for me\n",
" - I wanted poliastro to be super easy to install by avoiding the \"two language\" problem (this includes Windows)\n",
" - The native debugger is broken https://github.com/cython/cython/issues/1717\n",
" - I really don't want to worry about some gore details\n",
"\n",
"I don't have lots of experience with it, so I don't have solid arguments against it."
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### PyPy\n",
"\n",
"\n",
"\n",
"* PyPy is a super interesting alternative Python implementation https://pypy.org/\n",
"* I really really want to use it more, but there are some obstacles:\n",
" - The documentation is a bit poor, even the changelogs\n",
" - Lacks interest from the mainstream community (including snarky comments by Guido about \"nobody using it in production\")\n",
" - Difficult to install libraries with binary extensions, although conda-forge is helping! https://conda-forge.org/status/\n",
" - PyPy has several incompatibilities with manylinux1 wheels https://bitbucket.org/pypy/pypy/issues/2617/ although manylinux2010 is taking off!\n",
" - At the moment, requires rewriting NumPy-based code to use plain Python lists... impractical"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"### And a lot others...\n",
"\n",
"https://github.com/pfalcon/awesome-python-compilers/\n",
"\n",
"The number of projects keeps growing"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"## All of a sudden, in 2012...\n",
"\n",
"\n",
"\n",
"https://twitter.com/teoliphant/status/235789560678858752"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "slide"
}
},
"source": [
"# What is numba?\n",
"\n",
"\n",
"\n",
"> Numba is an open source JIT compiler that translates a subset of Python and NumPy code into fast machine code.\n",
"\n",
"* Latest stable version (at the time of writing) 0.48, (0.49.0rc1 tagged 13 days ago)\n",
"* Documentation https://numba.pydata.org/numba-doc/latest/index.html\n",
"* BSD-2 License\n",
"* Easy to install:\n",
"\n",
"```\n",
"$ pip install numba\n",
"$ conda install numba [--channel conda-forge]\n",
"```"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Example 1: Monte Carlo for $\\pi$\n",
"\n",
"\n",
"\n",
"(From https://commons.wikimedia.org/wiki/File:Pi_30K.gif)"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"import random\n",
"\n",
"def monte_carlo_pi(nsamples):\n",
" acc = 0\n",
" for i in range(nsamples):\n",
" x = random.random()\n",
" y = random.random()\n",
" if (x ** 2 + y ** 2) < 1.0:\n",
" acc += 1\n",
" return 4.0 * acc / nsamples"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"3.56"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"monte_carlo_pi(100)"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"data": {
"text/plain": [
"3.14984"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"monte_carlo_pi(100_000)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"5.38 s ± 708 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
]
}
],
"source": [
"%timeit monte_carlo_pi(10_000_000) # Slow!"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"from numba import jit"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"monte_carlo_pi_fast = jit(monte_carlo_pi)"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"272 ms ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
]
}
],
"source": [
"%timeit -n1 -r1 monte_carlo_pi_fast(100)"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"152 ms ± 8.73 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
]
}
],
"source": [
"%timeit monte_carlo_pi_fast(10_000_000) # 40x faster!"
]
},
{
"cell_type": "markdown",
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"source": [
"## Example 2: Plate deflection\n",
"\n",
"$$\n",
"w(\\xi, \\eta) = \\frac{4 P_c}{\\pi^4 D L_x L_y} \\sum_{m=1}^\\infty \\sum_{n=1}^\\infty \\frac{\\sin{\\frac{m \\pi \\xi}{L_x}} \\sin \\frac{n \\pi \\eta}{L_y}}{\\left(\\left(\\frac{m}{L_x}\\right)^2 + \\left(\\frac{n}{L_y}\\right)^2\\right)^2}\n",
"$$\n",
"\n",
"(From http://www.efunda.com/formulae/solid_mechanics/plates/calculators/SSSS_PPoint.cfm)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"import numpy as np\n",
"from numpy import pi, sin"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [],
"source": [
"@jit\n",
"def a_mn_point(P, a, b, xi, eta, mm, nn):\n",
" \"\"\"Navier series coefficient for concentrated load.\"\"\"\n",
" return 4 * P * sin(mm * pi * xi / a) * sin(nn * pi * eta / b) / (a * b)"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [],
"source": [
"@jit\n",
"def plate_displacement(xx, yy, ww, a, b, P, xi, eta, D, max_m, max_n):\n",
" max_i, max_j = ww.shape\n",
" for mm in range(1, max_m):\n",
" for nn in range(1, max_n):\n",
" for ii in range(max_i):\n",
" for jj in range(max_j):\n",
" a_mn = a_mn_point(P, a, b, xi, eta, mm, nn)\n",
" ww[ii, jj] += (\n",
" a_mn\n",
" / (mm ** 2 / a ** 2 + nn ** 2 / b ** 2) ** 2\n",
" * sin(mm * pi * xx[ii, jj] / a)\n",
" * sin(nn * pi * yy[ii, jj] / b)\n",
" / (pi ** 4 * D)\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"# Plate geometry\n",
"a = 1.0 # m\n",
"b = 1.0 # m\n",
"h = 50e-3 # m\n",
"\n",
"# Material properties\n",
"E = 69e9 # Pa\n",
"nu = 0.35\n",
"\n",
"# Series terms\n",
"max_m = 16\n",
"max_n = 16\n",
"\n",
"# Computation points\n",
"# NOTE: With an odd number of points the center of the place is included in\n",
"# the grid\n",
"NUM_POINTS = 101\n",
"\n",
"# Load\n",
"P = -10e3 # N\n",
"xi = a / 2\n",
"eta = a / 2\n",
"\n",
"# Flexural rigidity\n",
"D = h**3 * E / (12 * (1 - nu**2))\n",
"\n",
"# Set up domain\n",
"x = np.linspace(0, a, num=NUM_POINTS)\n",
"y = np.linspace(0, b, num=NUM_POINTS)\n",
"xx, yy = np.meshgrid(x, y)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"23.2 s ± 0 ns per loop (mean ± std. dev. of 1 run, 1 loop each)\n"
]
}
],
"source": [
"# Compute displacement field\n",
"ww = np.zeros_like(xx)\n",
"%timeit -n1 -r1 plate_displacement.py_func(xx, yy, ww, a, b, P, xi, eta, D, max_m, max_n)"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"slideshow": {
"slide_type": "fragment"
}
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"132 ms ± 3.11 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
]
}
],
"source": [
"%timeit plate_displacement(xx, yy, ww, a, b, P, xi, eta, D, max_m, max_n) # 100 times faster!"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"slideshow": {
"slide_type": "skip"
}
},
"outputs": [],
"source": [
"# https://github.com/plotly/plotly.py/issues/1664#issuecomment-511773518\n",
"import plotly.graph_objects as go\n",
"import plotly.io as pio\n",
"\n",
"# Set default renderer\n",
"pio.renderers.default = \"notebook_connected\"\n",
"\n",
"# Set default template\n",
"pio.templates[\"slides\"] = go.layout.Template(layout=dict(width=800, height=550))\n",
"pio.templates.default = \"plotly+slides\""
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"slideshow": {
"slide_type": "subslide"
}
},
"outputs": [
{
"data": {
"text/html": [
" \n",
" "
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
" \n",
" \n",
"
\n",
" \n",
"